#remotes::install_github("ropensci/rnoaa")
#remotes::install_github("ropensci/chirps")
#install.packages("GSODR")
#install.packages("AICcmodavg")The Effects of Climate Change on the Breeding Success and Timing of African White-backed Vultures
Load Libraries
library(tidyverse)── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 4.0.1 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.1
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(lubridate)
library(janitor)
Attaching package: 'janitor'
The following objects are masked from 'package:stats':
chisq.test, fisher.testlibrary(MASS)
Attaching package: 'MASS'
The following object is masked from 'package:dplyr':
selectlibrary(patchwork)
Attaching package: 'patchwork'
The following object is masked from 'package:MASS':
arealibrary(stringr)
library(tibble)
library(dplyr)
library(knitr)
library(ggplot2)
library(GGally)
library(tidyr)
library(chirps)
library(GSODR)
library(tibble)
library(knitr)
library(AICcmodavg)
library(lme4)Loading required package: Matrix
Attaching package: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
Attaching package: 'lme4'
The following object is masked from 'package:AICcmodavg':
checkConvOverdispersion Check
check_overdispersion <- function(model) {
rp <- residuals(model, type = "pearson")
disp <- sum(rp^2) / model$df.residual
disp
} vultures <- read.csv("~/Desktop/Vulture Project/Data/93_24 cleaned.csv")
# WeatherData <- read.csv("~/Desktop/Vulture Project/Data/Weather.csv")Weather Data
library(rnoaa)Registered S3 method overwritten by 'hoardr':
method from
print.cache_info httrThe rnoaa package will soon be retired and archived because the underlying APIs have changed dramatically. The package currently works but does not pull the most recent data in all cases. A noaaWeather package is planned as a replacement but the functions will not be interchangeable.library(tidyverse)
library(lubridate)
kim_station <- "684380-99999"
weather_raw <- get_GSOD(
station = kim_station,
years = 1993:2024
)WeatherData <- weather_raw |>
transmute(
YEAR,
MONTH,
DAY,
TEMP,
MAX,
MIN,
PRCP,
WDSP,
MXSPD,
I_HAIL
) |>
arrange(YEAR, MONTH, DAY)CHIRPS Rainfall for 2004
# coordinates (24.81 E, 28.62 S)
lonlat <- data.frame(
lon = 24.81,
lat = -28.62
)
kimberley_chirps_2004 <- get_chirps(
object = lonlat,
dates = c("2004-01-01", "2004-12-31"),
server = "ClimateSERV",
as.matrix = FALSE
)
Fetching data from ClimateSERV
Getting your request...glimpse(kimberley_chirps_2004)Rows: 366
Columns: 5
$ id <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ lon <dbl> 24.81, 24.81, 24.81, 24.81, 24.81, 24.81, 24.81, 24.81, 24.81, …
$ lat <dbl> -28.62, -28.62, -28.62, -28.62, -28.62, -28.62, -28.62, -28.62,…
$ date <date> 2004-01-01, 2004-01-02, 2004-01-03, 2004-01-04, 2004-01-05, 20…
$ chirps <dbl> 0.000000, 0.000000, 0.000000, 0.000000, 12.762877, 21.295383, 1…CHIRPS Rainfall for 2004 (Dronfield / Kimberley)
daily_chirps_2004 <- tibble(
date = as.Date(kimberley_chirps_2004$date),
rain_mm = as.numeric(kimberley_chirps_2004$chirps)
)
head(daily_chirps_2004) # A tibble: 6 × 2
date rain_mm
<date> <dbl>
1 2004-01-01 0
2 2004-01-02 0
3 2004-01-03 0
4 2004-01-04 0
5 2004-01-05 12.8
6 2004-01-06 21.3summary(daily_chirps_2004$rain_mm) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 0.00 0.00 1.26 0.00 28.76 annual_chirps_2004 <- daily_chirps_2004 |>
mutate(YEAR = year(date)) |>
group_by(YEAR) |>
summarise(
total_rain_2004 = sum(rain_mm, na.rm = TRUE),
max_rain_2004 = max(rain_mm, na.rm = TRUE),
rain_days_2004 = sum(rain_mm > 0, na.rm = TRUE),
.groups = "drop"
)
annual_chirps_2004# A tibble: 1 × 4
YEAR total_rain_2004 max_rain_2004 rain_days_2004
<dbl> <dbl> <dbl> <int>
1 2004 461. 28.8 54total_rain_2004 <- annual_chirps_2004$total_rain_2004
max_rain_2004 <- annual_chirps_2004$max_rain_2004
rain_days_2004 <- annual_chirps_2004$rain_days_2004
total_rain_2004[1] 461.1258max_rain_2004[1] 28.75726rain_days_2004[1] 54Weather Data
WeatherData_clean <- WeatherData |>
filter(
YEAR >= 1993,
YEAR <= 2024
) |>
mutate(
date = make_date(YEAR, MONTH, DAY)
) |>
# add CHIRPS daily rainfall into 2004
left_join(daily_chirps_2004, by = "date") |>
mutate(
PRCP = if_else(YEAR == 2004 & !is.na(rain_mm), rain_mm, PRCP)
) |>
dplyr::select(-rain_mm) vultures.clean <- vultures |>
clean_names() |> # ringing_date, laying_date
rename_with(~ gsub("_", ".", .x)) |> # ringing.date, laying.dater
mutate(
ringing.date = ringing.date |> as.character() |> str_trim(),
ringing.date = if_else(ringing.date %in% c("", "NA", "?"),
NA_character_, ringing.date),
ringing.date = ymd(gsub("/", "-", ringing.date)),
laying.date = laying.date |> as.character() |> str_trim(),
laying.date = if_else(laying.date %in% c("", "NA", "?"),
NA_character_, laying.date),
laying.date = ymd(gsub("/", "-", laying.date))
)Warning: There was 1 warning in `mutate()`.
ℹ In argument: `ringing.date = ymd(gsub("/", "-", ringing.date))`.
Caused by warning:
! 32 failed to parse.glimpse(vultures.clean)Rows: 2,657
Columns: 14
$ ringing.date <date> 1993-09-21, 1993-09-21, 1993-09-19, 1993-09-1…
$ tree.drn <chr> "1", "2", "3", "4", "5", "6", "8", "12", "13",…
$ southings <chr> " 28 38 43.8", "28 38 48.2", "28 38 55.6", "28…
$ eastings <chr> "24 47 14.5", "24 47 05.1", "24 47 15.4", "24 …
$ code <int> 1, 1, 1, 1, 1, 1, 1, 15, 1, 1, 1, 1, 3, 1, 5, …
$ ring <chr> "G19791", "G19792", "G19781", "G19780", "G1978…
$ colour.rings.r.top.down <chr> "M", "M", "M", "M", "M", "M", "M", "", "M", "M…
$ colour.rings.l.top.down <chr> "B, Y, Y", "B, Y, R", "B, W, G", "B, B, R", "B…
$ mass.g <int> 5900, 5000, 3900, 4950, 4950, 3800, 3550, NA, …
$ wing.length.mm <int> 475, 425, 260, 310, 290, 275, 210, NA, 245, 29…
$ age.at.ringing.days <int> 89, 80, 57, 64, 61, 59, 49, NA, 55, 61, 5, 92,…
$ value.at.ringing <int> 34233, 34233, 34231, 34231, 34231, 34231, 3423…
$ hatching.date <chr> "1993/06/24", "1993/07/03", "1993/07/24", "199…
$ laying.date <date> 1993-04-29, 1993-05-08, 1993-05-29, 1993-05-2…summary(vultures.clean) ringing.date tree.drn southings eastings
Min. :1993-09-19 Length:2657 Length:2657 Length:2657
1st Qu.:2004-06-06 Class :character Class :character Class :character
Median :2012-10-13 Mode :character Mode :character Mode :character
Mean :2011-08-31
3rd Qu.:2019-08-21
Max. :2024-10-18
NA's :32
code ring colour.rings.r.top.down
Min. : 1.000 Length:2657 Length:2657
1st Qu.: 1.000 Class :character Class :character
Median : 1.000 Mode :character Mode :character
Mean : 5.001
3rd Qu.: 6.000
Max. :138.000
NA's :4
colour.rings.l.top.down mass.g wing.length.mm age.at.ringing.days
Length:2657 Min. : 2 Min. : 25.0 Min. : 5.00
Class :character 1st Qu.:4200 1st Qu.:320.0 1st Qu.: 65.00
Mode :character Median :4800 Median :395.0 Median : 76.00
Mean :4610 Mean :375.7 Mean : 74.18
3rd Qu.:5200 3rd Qu.:460.0 3rd Qu.: 86.00
Max. :6700 Max. :570.0 Max. :300.00
NA's :1225 NA's :1192 NA's :1198
value.at.ringing hatching.date laying.date
Min. :34231 Length:2657 Min. :1900-01-02
1st Qu.:37541 Class :character 1st Qu.:2002-06-02
Median :40824 Mode :character Median :2011-05-13
Mean :40489 Mean :2010-05-21
3rd Qu.:43386 3rd Qu.:2018-05-26
Max. :45583 Max. :2024-08-11
NA's :1231 NA's :1197 colSums(is.na(vultures.clean)) ringing.date tree.drn southings
32 0 0
eastings code ring
0 4 0
colour.rings.r.top.down colour.rings.l.top.down mass.g
0 0 1225
wing.length.mm age.at.ringing.days value.at.ringing
1192 1198 1231
hatching.date laying.date
0 1197 Active Nests per year
active.nests.raw <- vultures.clean |>
mutate(year = year(ringing.date)) |>
group_by(year) |>
summarise(active_nests = n(), .groups = "drop")
ggplot(active.nests.raw, aes(x = year, y = active_nests)) +
geom_col(fill = "#89c5d6", alpha = 0.85) +
labs(x = "Year", y = "Number of active nests") +
theme_bw() +
theme(panel.grid = element_blank())Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_col()`).
Success proportion
failed.nests.raw <- vultures.clean |>
mutate(year = year(ringing.date)) |>
group_by(year) |>
summarise(
active_nests = n(),
success_nests = sum(!is.na(laying.date)),
success_prop = success_nests / active_nests,
.groups = "drop"
)
ggplot(failed.nests.raw, aes(x = success_prop)) +
geom_histogram(bins = 15, fill = "#89c5d6") +
labs(x = "Breeding success proportion", y = "Count") +
theme_bw() +
theme(panel.grid = element_blank())
Annual Weather Summaries
weather_yearly <- WeatherData_clean |>
group_by(YEAR) |>
summarise(
mean_temp = mean(TEMP, na.rm = TRUE),
max_temp = max(MAX, na.rm = TRUE),
min_temp = min(MIN, na.rm = TRUE),
total_rain = sum(PRCP, na.rm = TRUE),
max_rain = max(PRCP, na.rm = TRUE),
rain_days = sum(PRCP > 0, na.rm = TRUE),
mean_wind = mean(WDSP, na.rm = TRUE),
max_wind = max(MXSPD, na.rm = TRUE),
hail_days = sum(I_HAIL, na.rm = TRUE),
.groups = "drop"
)
# Check the corrected 2004 rainfall
weather_yearly[weather_yearly$YEAR == 2004,
c("YEAR", "total_rain", "max_rain", "rain_days")]# A tibble: 1 × 4
YEAR total_rain max_rain rain_days
<int> <dbl> <dbl> <int>
1 2004 461. 28.8 54Plot Annual Rainfall
ggplot(weather_yearly, aes(x = YEAR, y = total_rain)) +
geom_line() +
geom_point() +
theme_bw() +
labs(x = "Year", y = "Total annual rainfall (mm)")
Rain over time
library(grid)
rain_yearly <- weather_yearly
rain_long <- rain_yearly |>
pivot_longer(
cols = c(total_rain, max_rain, rain_days),
names_to = "metric",
values_to = "value"
) |>
mutate(
metric = factor(
metric,
levels = c("total_rain", "max_rain", "rain_days")
),
metric_label = case_when(
metric == "total_rain" ~ "Total annual rainfall (mm)",
metric == "max_rain" ~ "Max daily rainfall (mm)",
metric == "rain_days" ~ "Number of rainy days (> 0 mm)"
),
metric_label = factor(
metric_label,
levels = c(
"Total annual rainfall (mm)",
"Max daily rainfall (mm)",
"Number of rainy days (> 0 mm)"
)
)
)
ggplot(rain_long, aes(x = YEAR, y = value)) +
geom_line(linewidth = 0.7) +
geom_point(size = 1.3) +
facet_wrap(~ metric_label, ncol = 1, scales = "free_y") +
theme_classic() +
labs(
x = "Year",
y = NULL
) +
theme(
strip.text = element_text(size = 9, face = "bold"),
axis.title.x = element_text(size = 10),
axis.text = element_text(size = 8),
panel.spacing = unit(0.7, "lines")
)
Check min temp
# Yearly
ggplot(weather_yearly, aes(x = YEAR, y = min_temp)) +
geom_line() +
geom_point() +
theme_bw() +
labs(x = "Year", y = "Annual minimum temperature (°C)")
# Daily
weather_daily <- WeatherData_clean |>
group_by(YEAR, MONTH, DAY) |>
summarise(
mean_temp = mean(TEMP, na.rm = TRUE),
max_temp = max(MAX, na.rm = TRUE),
min_temp = min(MIN, na.rm = TRUE),
total_rain = sum(PRCP, na.rm = TRUE),
max_rain = max(PRCP, na.rm = TRUE),
rain_days = sum(PRCP > 0, na.rm = TRUE),
mean_wind = mean(WDSP, na.rm = TRUE),
max_wind = max(MXSPD, na.rm = TRUE),
hail_days = sum(I_HAIL, na.rm = TRUE),
.groups = "drop"
)Warning: There were 291 warnings in `summarise()`.
The first warning was:
ℹ In argument: `max_temp = max(MAX, na.rm = TRUE)`.
ℹ In group 9935: `YEAR = 2020`, `MONTH = 6`, `DAY = 15`.
Caused by warning in `max()`:
! no non-missing arguments to max; returning -Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 290 remaining warnings.temp_daily <- weather_daily |>
mutate(
min_temp = ifelse(is.infinite(min_temp), NA_real_, min_temp),
max_temp = ifelse(is.infinite(max_temp), NA_real_, max_temp)
)
summary(temp_daily$min_temp) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-9.900 4.500 10.750 9.818 15.500 30.000 9 summary(temp_daily$max_temp) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
5.50 23.70 28.60 28.18 32.90 45.70 1 Yearly temp summaries
weather_yearly_temp <- weather_daily |>
filter(YEAR >= 1993, YEAR <= 2024) |>
group_by(YEAR) |>
summarise(
mean_temp = mean(mean_temp, na.rm = TRUE),
max_temp = max(max_temp, na.rm = TRUE),
min_temp = min(min_temp, na.rm = TRUE),
.groups = "drop"
)
summary(weather_yearly_temp$min_temp) Min. 1st Qu. Median Mean 3rd Qu. Max.
-9.900 -7.425 -6.450 -6.356 -5.450 -3.300 summary(weather_yearly_temp$max_temp) Min. 1st Qu. Median Mean 3rd Qu. Max.
37.20 38.35 39.50 40.01 41.17 45.70 Temps over time
min_temp <- ggplot(weather_yearly_temp, aes(YEAR, min_temp)) +
geom_line() +
geom_point() +
theme_bw() +
labs(x = "Year", y = "Annual minimum temperature (°C)")
max_temp <- ggplot(weather_yearly_temp, aes(x = YEAR, y = max_temp)) +
geom_line() +
geom_point() +
theme_bw() +
labs(
x = "Year",
y = "Annual maximum temperature (°C)",
title = "Annual Maximum Temperature (1993–2024)"
)
mean_temp <- ggplot(weather_yearly_temp, aes(x = YEAR, y = mean_temp)) +
geom_line() +
geom_point() +
theme_bw() +
labs(
x = "Year",
y = "Annual mean temperature (°C)",
title = "Annual Mean Temperature (1993–2024)"
)temp_yearly_long <- weather_yearly_temp |>
pivot_longer(
cols = c(max_temp, min_temp, mean_temp),
names_to = "type",
values_to = "value"
) |>
mutate(
type = factor(
type,
levels = c("max_temp", "min_temp", "mean_temp"),
labels = c("Annual Maximum Temperature (°C)",
"Annual Minimum Temperature (°C)",
"Annual Mean Temperature (°C)")
)
)
ggplot(temp_yearly_long, aes(x = YEAR, y = value)) +
geom_line(size = 0.7) +
geom_point(size = 1.5) +
facet_wrap(~ type, scales = "free_y", ncol = 1) +
theme_classic(base_size = 11) +
theme(
strip.text = element_text(size = 10, face = "bold"),
axis.title = element_text(size = 9),
axis.text = element_text(size = 8),
panel.spacing = unit(0.8, "lines")
) +
labs(
x = "Year",
y = "",
title = "Annual Temperature Trends (1993–2024)"
)Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Lagged Weather summaries
weather_yearly_lag <- weather_yearly |>
arrange(YEAR) |>
mutate(
across(
c(mean_temp, max_temp, min_temp, total_rain, max_rain,
rain_days, mean_wind, max_wind, hail_days),
~ lag(.x, 1), .names = "lag_{.col}"
)
)Daily weather with lagged variables
weather_lag <- WeatherData_clean |>
mutate(
mean_temp = readr::parse_number(as.character(TEMP)),
max_temp = readr::parse_number(as.character(MAX)),
min_temp = readr::parse_number(as.character(MIN)),
mean_wind = readr::parse_number(as.character(WDSP)),
max_wind = readr::parse_number(as.character(MXSPD)),
rain = readr::parse_number(as.character(PRCP)),
rain_day = as.integer(replace_na(rain, 0) > 0),
hail = I_HAIL
) |>
arrange(date)Check for weather correlation
weather_yearly |>
dplyr::select(
mean_temp, max_temp, min_temp,
total_rain, max_rain, rain_days,
mean_wind, max_wind, hail_days
)# A tibble: 32 × 9
mean_temp max_temp min_temp total_rain max_rain rain_days mean_wind max_wind
<dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
1 18.4 39.1 -7.8 568. 111 67 4.05 15.4
2 17.8 39.6 -8 378. 37.3 39 3.75 13.9
3 18.5 39.1 -7.2 567. 81.3 70 3.79 18
4 17.7 37.5 -6.5 516. 39.1 79 3.98 14.9
5 18.5 40.9 -3.3 381. 33.0 65 3.76 13.9
6 18.7 37.8 -6 550. 32 84 4.00 20.6
7 19.4 37.9 -5 354. 42.2 62 3.95 28.3
8 18.0 37.2 -8.1 578. 78.5 88 3.61 20.6
9 19.0 39.9 -6.4 548. 27.4 105 3.78 27.8
10 19.4 43.6 -6.6 579. 44.7 83 3.71 24.2
# ℹ 22 more rows
# ℹ 1 more variable: hail_days <dbl>cor_mat <- weather_yearly |>
dplyr::select(mean_temp, max_temp, min_temp, total_rain, max_rain,
rain_days, mean_wind, max_wind, hail_days) |>
cor(use = "pairwise.complete.obs")
cor_df <- as.data.frame(cor_mat) |>
rownames_to_column("var1") |>
pivot_longer(-var1, names_to = "var2", values_to = "cor")
ggplot(cor_df, aes(x = var1, y = var2, fill = cor)) +
geom_tile() +
scale_fill_gradient2(
low = "#b2182b", mid = "white", high = "#2166ac",
midpoint = 0, limits = c(-1, 1)
) +
labs(x = "", y = "", fill = "Correlation") +
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
panel.grid = element_blank()
)
All weather across years
weather_yearly |>
pivot_longer(-YEAR) |>
ggplot(aes(x = YEAR, y = value)) +
geom_line(colour = "#89c5d6") +
facet_wrap(~ name, scales = "free_y") +
labs(x = "Year", y = "Value") +
theme_bw() +
theme(panel.grid = element_blank())
# Effort Data
active.nests <- vultures.clean |>
mutate(year = year(ringing.date)) |>
group_by(year) |>
summarise(active_nests = n(), .groups = "drop")
effort_weather <- active.nests |>
left_join(weather_yearly, by = c("year" = "YEAR"))Year only model
m_year <- glm.nb(active_nests ~ year, data = active.nests)
summary(m_year)
Call:
glm.nb(formula = active_nests ~ year, data = active.nests, init.theta = 283.4337529,
link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -53.188260 4.927882 -10.79 <2e-16 ***
year 0.028659 0.002451 11.69 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(283.4338) family taken to be 1)
Null deviance: 173.368 on 31 degrees of freedom
Residual deviance: 34.003 on 30 degrees of freedom
(1 observation deleted due to missingness)
AIC: 246.56
Number of Fisher Scoring iterations: 1
Theta: 283
Std. Err.: 329
2 x log-likelihood: -240.558 check_overdispersion(m_year)[1] 1.067024Autocorrelation in residuals
res_effort <- residuals(m_year, type = "pearson")
acf(res_effort,
main = "ACF of residuals: breeding effort (active nests)")
Rainfall
effort_weather <- effort_weather |>
mutate(
z_total_rain = scale(total_rain),
z_max_rain = scale(max_rain),
z_rain_days = scale(rain_days)
)
m_rain_only <- glm.nb(
active_nests ~ z_total_rain + z_max_rain + z_rain_days,
data = effort_weather
)
summary(m_rain_only)
Call:
glm.nb(formula = active_nests ~ z_total_rain + z_max_rain + z_rain_days,
data = effort_weather, init.theta = 18.0307885, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.39736 0.04606 95.461 <2e-16 ***
z_total_rain -0.18087 0.07301 -2.477 0.0132 *
z_max_rain -0.01308 0.05642 -0.232 0.8167
z_rain_days 0.11469 0.06637 1.728 0.0840 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(18.0308) family taken to be 1)
Null deviance: 41.511 on 31 degrees of freedom
Residual deviance: 32.342 on 28 degrees of freedom
(1 observation deleted due to missingness)
AIC: 294.97
Number of Fisher Scoring iterations: 1
Theta: 18.03
Std. Err.: 5.51
2 x log-likelihood: -284.973 Temp
effort_weather <- effort_weather |>
mutate(
z_mean_temp = scale(mean_temp),
z_max_temp = scale(max_temp),
z_min_temp = scale(min_temp)
)
max_temp_mod <- glm.nb(
active_nests ~ z_max_temp,
data = effort_weather
)
summary(max_temp_mod)
Call:
glm.nb(formula = active_nests ~ z_max_temp, data = effort_weather,
init.theta = 15.77748683, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.40154 0.04864 90.494 <2e-16 ***
z_max_temp 0.10583 0.04911 2.155 0.0312 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(15.7775) family taken to be 1)
Null deviance: 37.216 on 31 degrees of freedom
Residual deviance: 32.494 on 30 degrees of freedom
(1 observation deleted due to missingness)
AIC: 294.66
Number of Fisher Scoring iterations: 1
Theta: 15.78
Std. Err.: 4.72
2 x log-likelihood: -288.663 min_temp_mod <- glm.nb(
active_nests ~ z_min_temp,
data = effort_weather
)
summary(min_temp_mod)
Call:
glm.nb(formula = active_nests ~ z_min_temp, data = effort_weather,
init.theta = 13.63264473, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.40648 0.05171 85.218 <2e-16 ***
z_min_temp 0.03703 0.05254 0.705 0.481
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(13.6326) family taken to be 1)
Null deviance: 32.929 on 31 degrees of freedom
Residual deviance: 32.466 on 30 degrees of freedom
(1 observation deleted due to missingness)
AIC: 298.61
Number of Fisher Scoring iterations: 1
Theta: 13.63
Std. Err.: 3.97
2 x log-likelihood: -292.606 mean_temp_mod <- glm.nb(
active_nests ~ z_mean_temp,
data = effort_weather
)
summary(mean_temp_mod)
Call:
glm.nb(formula = active_nests ~ z_mean_temp, data = effort_weather,
init.theta = 13.63790192, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.40648 0.05170 85.233 <2e-16 ***
z_mean_temp 0.03478 0.05251 0.662 0.508
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(13.6379) family taken to be 1)
Null deviance: 32.939 on 31 degrees of freedom
Residual deviance: 32.476 on 30 degrees of freedom
(1 observation deleted due to missingness)
AIC: 298.61
Number of Fisher Scoring iterations: 1
Theta: 13.64
Std. Err.: 3.97
2 x log-likelihood: -292.605 Wind
effort_weather <- effort_weather |>
mutate(
z_mean_wind = scale(mean_wind),
z_max_wind = scale(max_wind)
)
m_wind_only <- glm.nb(
active_nests ~ z_mean_wind + z_max_wind,
data = effort_weather
)
summary(m_wind_only)
Call:
glm.nb(formula = active_nests ~ z_mean_wind + z_max_wind, data = effort_weather,
init.theta = 15.33069655, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.40243 0.04922 89.438 <2e-16 ***
z_mean_wind -0.06706 0.05215 -1.286 0.199
z_max_wind -0.05428 0.05231 -1.038 0.299
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(15.3307) family taken to be 1)
Null deviance: 36.340 on 31 degrees of freedom
Residual deviance: 32.489 on 29 degrees of freedom
(1 observation deleted due to missingness)
AIC: 297.43
Number of Fisher Scoring iterations: 1
Theta: 15.33
Std. Err.: 4.56
2 x log-likelihood: -289.431 Hail
effort_weather <- effort_weather |>
mutate(z_hail_days = scale(hail_days))
m_hail_only <- glm.nb(
active_nests ~ z_hail_days,
data = effort_weather
)
summary(m_hail_only)
Call:
glm.nb(formula = active_nests ~ z_hail_days, data = effort_weather,
init.theta = 15.4472834, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.40226 0.04907 89.721 <2e-16 ***
z_hail_days 0.09695 0.04926 1.968 0.049 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(15.4473) family taken to be 1)
Null deviance: 36.569 on 31 degrees of freedom
Residual deviance: 32.511 on 30 degrees of freedom
(1 observation deleted due to missingness)
AIC: 295.25
Number of Fisher Scoring iterations: 1
Theta: 15.45
Std. Err.: 4.60
2 x log-likelihood: -289.249 Active nests across years
active.nests <- vultures.clean |>
mutate(year = year(ringing.date)) |>
group_by(year) |>
summarise(active_nests = n(), .groups = "drop")
m_year <- glm.nb(active_nests ~ year, data = active.nests)
active.nests_plot <- active.nests |>
mutate(pred_year = predict(m_year, newdata = active.nests, type = "response"))
ggplot(active.nests_plot, aes(x = year, y = active_nests)) +
geom_col(fill = "#89c5d6", alpha = 0.85) +
geom_line(aes(y = pred_year), linewidth = 1.1, colour = "black") +
labs(
x = "Year",
y = "Number of active nests"
) +
theme_bw() +
theme(panel.grid = element_blank())Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_col()`).Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_line()`).
Effects of max temp on effort
effort_weather <- active.nests |>
left_join(weather_yearly, by = c("year" = "YEAR")) |>
mutate(
z_total_rain = as.numeric(scale(total_rain)),
z_max_rain = as.numeric(scale(max_rain)),
z_rain_days = as.numeric(scale(rain_days)),
z_mean_temp = as.numeric(scale(mean_temp)),
z_max_temp = as.numeric(scale(max_temp)),
z_min_temp = as.numeric(scale(min_temp)),
z_mean_wind = as.numeric(scale(mean_wind)),
z_max_wind = as.numeric(scale(max_wind)),
z_hail_days = as.numeric(scale(hail_days))
)
m_eff_maxtemp <- glm.nb(
active_nests ~ z_max_temp,
data = effort_weather
)
summary(m_eff_maxtemp)
Call:
glm.nb(formula = active_nests ~ z_max_temp, data = effort_weather,
init.theta = 15.77748683, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.40154 0.04864 90.494 <2e-16 ***
z_max_temp 0.10583 0.04911 2.155 0.0312 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(15.7775) family taken to be 1)
Null deviance: 37.216 on 31 degrees of freedom
Residual deviance: 32.494 on 30 degrees of freedom
(1 observation deleted due to missingness)
AIC: 294.66
Number of Fisher Scoring iterations: 1
Theta: 15.78
Std. Err.: 4.72
2 x log-likelihood: -288.663 new_eff_temp <- tibble(
z_max_temp = seq(
min(effort_weather$z_max_temp, na.rm = TRUE),
max(effort_weather$z_max_temp, na.rm = TRUE),
length.out = 100
)
) |>
mutate(
pred_nests = predict(
m_eff_maxtemp,
newdata = cur_data(),
type = "response"
)
)Warning: There was 1 warning in `mutate()`.
ℹ In argument: `pred_nests = predict(m_eff_maxtemp, newdata = cur_data(), type
= "response")`.
Caused by warning:
! `cur_data()` was deprecated in dplyr 1.1.0.
ℹ Please use `pick()` instead.ggplot(effort_weather, aes(x = z_max_temp, y = active_nests)) +
geom_point(alpha = 0.7) +
geom_line(
data = new_eff_temp,
aes(x = z_max_temp, y = pred_nests),
linewidth = 1.1
) +
labs(
x = "Annual maximum temperature (scaled)",
y = "Number of active nests"
) +
theme_bw() +
theme(panel.grid = element_blank())Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_point()`).
acf(active.nests$active_nests)
Success DF
failed.nests <- vultures.clean |>
mutate(year = year(ringing.date)) |>
filter(!is.na(year)) |>
group_by(year) |>
summarise(
active_nests = n(),
success_nests = sum(!is.na(laying.date)),
failed_nests = active_nests - success_nests,
success_prop = success_nests / active_nests,
.groups = "drop"
) |>
filter(active_nests > 0)Year only model
success_year <- glm(
cbind(success_nests, failed_nests) ~ year,
family = binomial,
data = failed.nests
)
summary(success_year)
Call:
glm(formula = cbind(success_nests, failed_nests) ~ year, family = binomial,
data = failed.nests)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 45.434641 8.858715 5.129 2.92e-07 ***
year -0.022505 0.004405 -5.109 3.24e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 107.689 on 31 degrees of freedom
Residual deviance: 81.283 on 30 degrees of freedom
AIC: 237.58
Number of Fisher Scoring iterations: 3check_overdispersion(success_year)[1] 2.670817Rain
success_df <- failed.nests |>
left_join(weather_yearly, by = c("year" = "YEAR")) |>
mutate(
z_total_rain = as.numeric(scale(total_rain)),
z_max_rain = as.numeric(scale(max_rain)),
z_rain_days = as.numeric(scale(rain_days))
)
success_rain <- glm(
cbind(success_nests, failed_nests) ~
z_total_rain + z_max_rain + z_rain_days,
family = binomial,
data = success_df
)
summary(success_rain)
Call:
glm(formula = cbind(success_nests, failed_nests) ~ z_total_rain +
z_max_rain + z_rain_days, family = binomial, data = success_df)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.18710 0.03966 4.717 2.39e-06 ***
z_total_rain 0.08304 0.06426 1.292 0.1963
z_max_rain 0.02901 0.05086 0.570 0.5685
z_rain_days -0.09864 0.05967 -1.653 0.0983 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 107.69 on 31 degrees of freedom
Residual deviance: 102.75 on 28 degrees of freedom
AIC: 263.05
Number of Fisher Scoring iterations: 3check_overdispersion(success_rain)[1] 3.581363Temp
success_df <- success_df |>
mutate(
z_mean_temp = as.numeric(scale(mean_temp)),
z_max_temp = as.numeric(scale(max_temp)),
z_min_temp = as.numeric(scale(min_temp))
)
success_minT <- glm(
cbind(success_nests, failed_nests) ~ z_min_temp,
family = binomial,
data = success_df
)
summary(success_minT)
Call:
glm(formula = cbind(success_nests, failed_nests) ~ z_min_temp,
family = binomial, data = success_df)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.17203 0.03931 4.376 1.21e-05 ***
z_min_temp 0.16071 0.04191 3.834 0.000126 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 107.689 on 31 degrees of freedom
Residual deviance: 92.877 on 30 degrees of freedom
AIC: 249.18
Number of Fisher Scoring iterations: 3check_overdispersion(success_minT)[1] 3.031976success_maxT <- glm(
cbind(success_nests, failed_nests) ~ z_max_temp,
family = binomial,
data = success_df
)
summary(success_maxT)
Call:
glm(formula = cbind(success_nests, failed_nests) ~ z_max_temp,
family = binomial, data = success_df)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.18310 0.03944 4.643 3.44e-06 ***
z_max_temp -0.06103 0.03839 -1.590 0.112
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 107.69 on 31 degrees of freedom
Residual deviance: 105.16 on 30 degrees of freedom
AIC: 261.46
Number of Fisher Scoring iterations: 3check_overdispersion(success_maxT)[1] 3.433179success_meanT <- glm(
cbind(success_nests, failed_nests) ~ z_mean_temp,
family = binomial,
data = success_df
)
summary(success_meanT)
Call:
glm(formula = cbind(success_nests, failed_nests) ~ z_mean_temp,
family = binomial, data = success_df)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.17461 0.03922 4.452 8.51e-06 ***
z_mean_temp 0.05587 0.03856 1.449 0.147
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 107.69 on 31 degrees of freedom
Residual deviance: 105.59 on 30 degrees of freedom
AIC: 261.89
Number of Fisher Scoring iterations: 3check_overdispersion(success_meanT)[1] 3.4372WInd
success_df <- success_df |>
mutate(
z_mean_wind = as.numeric(scale(mean_wind)),
z_max_wind = as.numeric(scale(max_wind))
)
m_success_wind <- glm(
cbind(success_nests, failed_nests) ~ z_mean_wind + z_max_wind,
family = binomial,
data = success_df
)
summary(m_success_wind)
Call:
glm(formula = cbind(success_nests, failed_nests) ~ z_mean_wind +
z_max_wind, family = binomial, data = success_df)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.17646 0.03937 4.482 7.4e-06 ***
z_mean_wind 0.01436 0.04176 0.344 0.731
z_max_wind -0.01688 0.04351 -0.388 0.698
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 107.69 on 31 degrees of freedom
Residual deviance: 107.48 on 29 degrees of freedom
AIC: 265.78
Number of Fisher Scoring iterations: 3check_overdispersion(m_success_wind)[1] 3.625179Hail
success_df <- success_df |>
mutate(
z_hail_days = as.numeric(scale(hail_days))
)
m_success_hail <- glm(
cbind(success_nests, failed_nests) ~ z_hail_days,
family = binomial,
data = success_df
)
summary(m_success_hail)
Call:
glm(formula = cbind(success_nests, failed_nests) ~ z_hail_days,
family = binomial, data = success_df)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.18275 0.03940 4.638 3.52e-06 ***
z_hail_days -0.05984 0.03647 -1.641 0.101
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 107.69 on 31 degrees of freedom
Residual deviance: 105.00 on 30 degrees of freedom
AIC: 261.3
Number of Fisher Scoring iterations: 3check_overdispersion(m_success_hail)[1] 3.42445160 and 30 day pre laying period
timing_df <- vultures.clean |>
filter(!is.na(laying.date)) |>
mutate(
year = year(laying.date),
lay_DOY = yday(laying.date)
)
summarise_60 <- function(lay_date, weather) {
w <- weather |>
filter(date > (lay_date - days(60)), date <= lay_date)
tibble(
mean_temp_60 = mean(w$mean_temp, na.rm = TRUE),
max_temp_60 = mean(w$max_temp, na.rm = TRUE),
min_temp_60 = mean(w$min_temp, na.rm = TRUE),
total_rain_60 = sum(w$rain, na.rm = TRUE),
rain_days_60 = sum(w$rain_day, na.rm = TRUE),
mean_wind_60 = mean(w$mean_wind, na.rm = TRUE),
max_wind_60 = mean(w$max_wind, na.rm = TRUE),
hail_days_60 = sum(w$hail, na.rm = TRUE)
)
}
timing_weather_60 <- timing_df |>
mutate(win60 = purrr::map(laying.date, summarise_60, weather = weather_lag)) |>
unnest(win60)
summarise_30 <- function(lay_date, weather) {
w <- weather |>
filter(date > (lay_date - days(30)), date <= lay_date)
tibble(
mean_temp_30 = mean(w$mean_temp, na.rm = TRUE),
max_temp_30 = mean(w$max_temp, na.rm = TRUE),
min_temp_30 = mean(w$min_temp, na.rm = TRUE),
total_rain_30 = sum(w$rain, na.rm = TRUE),
rain_days_30 = sum(w$rain_day, na.rm = TRUE),
mean_wind_30 = mean(w$mean_wind, na.rm = TRUE),
max_wind_30 = mean(w$max_wind, na.rm = TRUE),
hail_days_30 = sum(w$hail, na.rm = TRUE)
)
}
timing_weather_30 <- timing_df |>
mutate(win30 = purrr::map(laying.date, summarise_30, weather = weather_lag)) |>
tidyr::unnest(win30)Temp
m_timing_temp_60 <- lm(
lay_DOY ~ mean_temp_60 + max_temp_60 + min_temp_60,
data = timing_weather_60
)
summary(m_timing_temp_60)
Call:
lm(formula = lay_DOY ~ mean_temp_60 + max_temp_60 + min_temp_60,
data = timing_weather_60)
Residuals:
Min 1Q Median 3Q Max
-88.511 -5.887 -1.339 4.323 145.827
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 230.2640 4.2240 54.514 < 2e-16 ***
mean_temp_60 2.6343 0.6348 4.150 3.52e-05 ***
max_temp_60 -3.0285 0.3587 -8.443 < 2e-16 ***
min_temp_60 -5.5674 0.3681 -15.123 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.13 on 1455 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.6752, Adjusted R-squared: 0.6745
F-statistic: 1008 on 3 and 1455 DF, p-value: < 2.2e-16Rain
m_timing_rain_60 <- lm(
lay_DOY ~ total_rain_60 + rain_days_60,
data = timing_weather_60
)
summary(m_timing_rain_60)
Call:
lm(formula = lay_DOY ~ total_rain_60 + rain_days_60, data = timing_weather_60)
Residuals:
Min 1Q Median 3Q Max
-166.729 -9.039 -2.407 7.930 134.078
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 168.72873 0.85555 197.217 < 2e-16 ***
total_rain_60 -0.06469 0.01242 -5.207 2.19e-07 ***
rain_days_60 -0.96694 0.09618 -10.054 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.69 on 1457 degrees of freedom
Multiple R-squared: 0.2568, Adjusted R-squared: 0.2557
F-statistic: 251.7 on 2 and 1457 DF, p-value: < 2.2e-16Hail
m_timing_rain_60 <- lm(
lay_DOY ~ total_rain_60 + rain_days_60,
data = timing_weather_60
)
summary(m_timing_rain_60)
Call:
lm(formula = lay_DOY ~ total_rain_60 + rain_days_60, data = timing_weather_60)
Residuals:
Min 1Q Median 3Q Max
-166.729 -9.039 -2.407 7.930 134.078
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 168.72873 0.85555 197.217 < 2e-16 ***
total_rain_60 -0.06469 0.01242 -5.207 2.19e-07 ***
rain_days_60 -0.96694 0.09618 -10.054 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.69 on 1457 degrees of freedom
Multiple R-squared: 0.2568, Adjusted R-squared: 0.2557
F-statistic: 251.7 on 2 and 1457 DF, p-value: < 2.2e-16Test
plot(failed.nests$year, failed.nests$success_prop)
Objective 1
objective1_table <- tribble(
~`Hypothesis / justification`, ~Model,
# ---- Temperature ----
"**Temp**", NA,
"Warmer years increase energetic cost of adult birds, reducing likelihood to breed.",
"active_nests ~ mean_temp",
"Extreme heat events cause thermal stress, reducing breeding effort (active cooling is energetically costly and may reduce adult condition).",
"active_nests ~ max_temp",
"Cold extremes increase thermoregulatory demand (often increasing time spent foraging), reducing effort put into breeding",
"active_nests ~ min_temp",
# ---- Rain ----
"**Rain**", NA,
"High rainfall reduces foraging efficiency and adult condition, reducing breeding effort.",
"active_nests ~ total_rain",
"Prolonged rainfall conditions reduce breeding effort by limiting foraging opportunity (fewer or weaker thermals for soaring)",
"active_nests ~ rain_days",
"Intense rainfall events decrease foraging efficiency and nest attendance, reducing breeding effort.",
"active_nests ~ max_rain",
"Previous year’s rainfall influences food availability and adult condition, affecting breeding effort in the following year.",
"active_nests ~ lag_total_rain",
# ---- Wind ----
"**Wind**", NA,
"Persistent windy conditions increase the cost of flying, limiting foraging and the energy available for reproduction, reducing breeding effort",
"active_nests ~ mean_wind",
"Extreme winds increase nest damage risk, reducing the likelihood that adults will attempt breeding (WBV often reuse nests from previous years) ",
"active_nests ~ max_wind",
# ---- Hail ----
"**Hail**", NA,
"Severe storm events increase disturbance and nest damage risk, reducing breeding effort",
"active_nests ~ hail_days",
# ---- Joint Effects ----
"**Joint effects**", NA,
"Cold stress combined with prolonged wet conditions jointly limit breeding effort",
"active_nests ~ min_temp + total_rain",
"Cold stress combined with intense rainfall events reduces breeding effort",
"active_nests ~ min_temp + max_rain",
"Cold extremes combined with extreme wind increase energetic costs and breeding risk, reducing breeding effort",
"active_nests ~ min_temp + max_wind",
"Severe storm exposure combined with extreme wind increases nest disturbance, reducing breeding effort",
"active_nests ~ hail_days + max_wind",
"Reduced thermals for soaring combined with high flight costs limit breeding effort",
"active_nests ~ max_rain + max_wind",
"Heat stress combined with prolonged wet conditions jointly reduce breeding effort by increasing energetic costs while limiting foraging",
"active_nests ~ max_temp + rain_days",
" High annual rainfall and frequent hail events each impose energetic and disturbance costs that reduce breeding effort.",
"active_nests ~ hail_days + rain_days",
# ---- Interactive Effects ----
"**Interactive effects**", NA,
"The negative effect of cold extremes on breeding effort is increased in persistently wet years.",
"active_nests ~ min_temp * rain_days",
"Cold-related energetic stress is amplified in windy years due to increased heat loss and flight costs",
"active_nests ~ min_temp * mean_wind",
"Adult condition resulting from the previous year interacts with current-year heat stress to influence breeding effort",
"active_nests ~ lag_total_rain * max_temp"
)
objective1_table_clean <- objective1_table |>
dplyr::mutate(Model = ifelse(is.na(Model), "", Model))
kable(
objective1_table_clean,
align = c("l", "l"),
caption = "Objective 1 candidate model set for breeding effort (active nests)"
)| Hypothesis / justification | Model |
|---|---|
| Temp | |
| Warmer years increase energetic cost of adult birds, reducing likelihood to breed. | active_nests ~ mean_temp |
| Extreme heat events cause thermal stress, reducing breeding effort (active cooling is energetically costly and may reduce adult condition). | active_nests ~ max_temp |
| Cold extremes increase thermoregulatory demand (often increasing time spent foraging), reducing effort put into breeding | active_nests ~ min_temp |
| Rain | |
| High rainfall reduces foraging efficiency and adult condition, reducing breeding effort. | active_nests ~ total_rain |
| Prolonged rainfall conditions reduce breeding effort by limiting foraging opportunity (fewer or weaker thermals for soaring) | active_nests ~ rain_days |
| Intense rainfall events decrease foraging efficiency and nest attendance, reducing breeding effort. | active_nests ~ max_rain |
| Previous year’s rainfall influences food availability and adult condition, affecting breeding effort in the following year. | active_nests ~ lag_total_rain |
| Wind | |
| Persistent windy conditions increase the cost of flying, limiting foraging and the energy available for reproduction, reducing breeding effort | active_nests ~ mean_wind |
| Extreme winds increase nest damage risk, reducing the likelihood that adults will attempt breeding (WBV often reuse nests from previous years) | active_nests ~ max_wind |
| Hail | |
| Severe storm events increase disturbance and nest damage risk, reducing breeding effort | active_nests ~ hail_days |
| Joint effects | |
| Cold stress combined with prolonged wet conditions jointly limit breeding effort | active_nests ~ min_temp + total_rain |
| Cold stress combined with intense rainfall events reduces breeding effort | active_nests ~ min_temp + max_rain |
| Cold extremes combined with extreme wind increase energetic costs and breeding risk, reducing breeding effort | active_nests ~ min_temp + max_wind |
| Severe storm exposure combined with extreme wind increases nest disturbance, reducing breeding effort | active_nests ~ hail_days + max_wind |
| Reduced thermals for soaring combined with high flight costs limit breeding effort | active_nests ~ max_rain + max_wind |
| Heat stress combined with prolonged wet conditions jointly reduce breeding effort by increasing energetic costs while limiting foraging | active_nests ~ max_temp + rain_days |
| High annual rainfall and frequent hail events each impose energetic and disturbance costs that reduce breeding effort. | active_nests ~ hail_days + rain_days |
| Interactive effects | |
| The negative effect of cold extremes on breeding effort is increased in persistently wet years. | active_nests ~ min_temp * rain_days |
| Cold-related energetic stress is amplified in windy years due to increased heat loss and flight costs | active_nests ~ min_temp * mean_wind |
| Adult condition resulting from the previous year interacts with current-year heat stress to influence breeding effort | active_nests ~ lag_total_rain * max_temp |
# Annual effort: number of active nests per year
active.nests <- vultures.clean |>
mutate(year = year(ringing.date)) |>
group_by(year) |>
summarise(active_nests = n(), .groups = "drop")
# Join to weather, create log effort response
effort_weather <- active.nests |>
left_join(weather_yearly_lag, by = c("year" = "YEAR")) |>
filter(!is.na(lag_total_rain)) |>
mutate(
log_active_nests = log(active_nests)
)# checking residuals
par(mfrow = c(1, 2))Fitingt then no-year models
# Baseline
m_null <- lm(log_active_nests ~ 1, data = effort_weather)
# Single predictor models
m_mean_temp <- lm(log_active_nests ~ mean_temp, data = effort_weather)
m_max_temp <- lm(log_active_nests ~ max_temp, data = effort_weather)
m_min_temp <- lm(log_active_nests ~ min_temp, data = effort_weather)
m_total_rain <- lm(log_active_nests ~ total_rain, data = effort_weather)
m_rain_days <- lm(log_active_nests ~ rain_days, data = effort_weather)
m_max_rain <- lm(log_active_nests ~ max_rain, data = effort_weather)
m_lag_total_rain <- lm(log_active_nests ~ lag_total_rain, data = effort_weather)
m_mean_wind <- lm(log_active_nests ~ mean_wind, data = effort_weather)
m_max_wind <- lm(log_active_nests ~ max_wind, data = effort_weather)
m_hail_days <- lm(log_active_nests ~ hail_days, data = effort_weather)
# Joint effects
m_minT_totalR <- lm(log_active_nests ~ min_temp + total_rain, data = effort_weather)
m_minT_maxR <- lm(log_active_nests ~ min_temp + max_rain, data = effort_weather)
m_minT_maxW <- lm(log_active_nests ~ min_temp + max_wind, data = effort_weather)
m_hail_maxW <- lm(log_active_nests ~ hail_days + max_wind, data = effort_weather)
m_maxR_maxW <- lm(log_active_nests ~ max_rain + max_wind, data = effort_weather)
m_maxT_rainDays <- lm(log_active_nests ~ max_temp + rain_days, data = effort_weather)
m_hail_rainDays <- lm(log_active_nests ~ hail_days + rain_days, data = effort_weather)
# Interactions
m_minT_x_rainDays <- lm(log_active_nests ~ min_temp * rain_days, data = effort_weather)
m_minT_x_meanWind <- lm(log_active_nests ~ min_temp * mean_wind, data = effort_weather)
m_lagR_x_maxT <- lm(log_active_nests ~ lag_total_rain * max_temp, data = effort_weather)Fit with year models
# Baseline
m_year <- lm(log_active_nests ~ year, data = effort_weather)
# Single predictor models
m_mean_temp_year <- lm(log_active_nests ~ year + mean_temp, data = effort_weather)
m_max_temp_year <- lm(log_active_nests ~ year + max_temp, data = effort_weather)
m_min_temp_year <- lm(log_active_nests ~ year + min_temp, data = effort_weather)
m_total_rain_year <- lm(log_active_nests ~ year + total_rain, data = effort_weather)
m_rain_days_year <- lm(log_active_nests ~ year + rain_days, data = effort_weather)
m_max_rain_year <- lm(log_active_nests ~ year + max_rain, data = effort_weather)
m_lag_total_rain_year <- lm(log_active_nests ~ year + lag_total_rain, data = effort_weather)
m_mean_wind_year <- lm(log_active_nests ~ year + mean_wind, data = effort_weather)
m_max_wind_year <- lm(log_active_nests ~ year + max_wind, data = effort_weather)
m_hail_days_year <- lm(log_active_nests ~ year + hail_days, data = effort_weather)
# Joint effects
m_minT_totalR_year <- lm(log_active_nests ~ year + min_temp + total_rain, data = effort_weather)
m_minT_maxR_year <- lm(log_active_nests ~ year + min_temp + max_rain, data = effort_weather)
m_minT_maxW_year <- lm(log_active_nests ~ year + min_temp + max_wind, data = effort_weather)
m_hail_maxW_year <- lm(log_active_nests ~ year + hail_days + max_wind, data = effort_weather)
m_maxR_maxW_year <- lm(log_active_nests ~ year + max_rain + max_wind, data = effort_weather)
m_maxT_rainDays_year <- lm(log_active_nests ~ year + max_temp + rain_days, data = effort_weather)
m_hail_rainDays_year <- lm(log_active_nests ~ year + hail_days + rain_days, data = effort_weather)
# Interactions
m_minT_x_rainDays_year <- lm(log_active_nests ~ year + min_temp * rain_days, data = effort_weather)
m_minT_x_meanWind_year <- lm(log_active_nests ~ year + min_temp * mean_wind, data = effort_weather)
m_lagR_x_maxT_year <- lm(log_active_nests ~ year + lag_total_rain * max_temp, data = effort_weather)Two separate AICc model selection runs
# AIC 1: without year
cand_effort_no_year <- list(
null = m_null,
mean_temp = m_mean_temp,
max_temp = m_max_temp,
min_temp = m_min_temp,
total_rain = m_total_rain,
rain_days = m_rain_days,
max_rain = m_max_rain,
lag_total_rain = m_lag_total_rain,
mean_wind = m_mean_wind,
max_wind = m_max_wind,
hail_days = m_hail_days,
minT_totalR = m_minT_totalR,
minT_maxR = m_minT_maxR,
minT_maxW = m_minT_maxW,
hail_maxW = m_hail_maxW,
maxR_maxW = m_maxR_maxW,
maxT_rainDays = m_maxT_rainDays,
hail_rainDays = m_hail_rainDays,
minT_x_rainDays = m_minT_x_rainDays,
minT_x_meanWind = m_minT_x_meanWind,
lagR_x_maxT = m_lagR_x_maxT
)
aic_effort_no_year <- aictab(cand.set = cand_effort_no_year, modnames = names(cand_effort_no_year))
aic_effort_no_year
Model selection based on AICc:
K AICc Delta_AICc AICcWt Cum.Wt LL
hail_days 3 12.04 0.00 0.17 0.17 -2.57
max_temp 3 12.77 0.74 0.12 0.29 -2.94
total_rain 3 12.83 0.79 0.12 0.41 -2.97
hail_maxW 4 13.16 1.12 0.10 0.51 -1.81
null 2 14.26 2.22 0.06 0.57 -4.92
maxT_rainDays 4 14.37 2.33 0.05 0.62 -2.41
max_wind 3 14.39 2.35 0.05 0.68 -3.75
hail_rainDays 4 14.44 2.40 0.05 0.73 -2.45
mean_wind 3 14.84 2.81 0.04 0.77 -3.98
lag_total_rain 3 15.07 3.03 0.04 0.81 -4.09
minT_totalR 4 15.45 3.41 0.03 0.84 -2.96
lagR_x_maxT 5 15.71 3.67 0.03 0.87 -1.65
max_rain 3 15.74 3.70 0.03 0.89 -4.43
maxR_maxW 4 15.84 3.81 0.03 0.92 -3.15
mean_temp 3 16.46 4.42 0.02 0.94 -4.79
min_temp 3 16.61 4.58 0.02 0.96 -4.86
rain_days 3 16.72 4.68 0.02 0.97 -4.92
minT_maxW 4 17.03 5.00 0.01 0.99 -3.75
minT_maxR 4 18.24 6.20 0.01 1.00 -4.35
minT_x_meanWind 5 20.27 8.23 0.00 1.00 -3.93
minT_x_rainDays 5 21.95 9.92 0.00 1.00 -4.78# AIC 2: with year
cand_effort_with_year <- list(
year = m_year,
mean_temp = m_mean_temp_year,
max_temp = m_max_temp_year,
min_temp = m_min_temp_year,
total_rain = m_total_rain_year,
rain_days = m_rain_days_year,
max_rain = m_max_rain_year,
lag_total_rain = m_lag_total_rain_year,
mean_wind = m_mean_wind_year,
max_wind = m_max_wind_year,
hail_days = m_hail_days_year,
minT_totalR = m_minT_totalR_year,
minT_maxR = m_minT_maxR_year,
minT_maxW = m_minT_maxW_year,
hail_maxW = m_hail_maxW_year,
maxR_maxW = m_maxR_maxW_year,
maxT_rainDays = m_maxT_rainDays_year,
hail_rainDays = m_hail_rainDays_year,
minT_x_rainDays = m_minT_x_rainDays_year,
minT_x_meanWind = m_minT_x_meanWind_year,
lagR_x_maxT = m_lagR_x_maxT_year
)
aic_effort_with_year <- aictab(cand.set = cand_effort_with_year, modnames = names(cand_effort_with_year))
aic_effort_with_year
Model selection based on AICc:
K AICc Delta_AICc AICcWt Cum.Wt LL
max_rain 4 -29.07 0.00 0.19 0.19 19.31
year 3 -28.43 0.64 0.14 0.33 17.66
total_rain 4 -27.49 1.58 0.09 0.42 18.51
minT_maxR 5 -27.44 1.63 0.09 0.51 19.92
min_temp 4 -26.68 2.39 0.06 0.56 18.11
mean_temp 4 -26.59 2.48 0.06 0.62 18.06
max_temp 4 -26.42 2.65 0.05 0.67 17.98
maxR_maxW 5 -26.41 2.66 0.05 0.72 19.40
max_wind 4 -25.90 3.18 0.04 0.76 17.72
hail_days 4 -25.88 3.19 0.04 0.80 17.71
mean_wind 4 -25.87 3.21 0.04 0.84 17.70
lag_total_rain 4 -25.81 3.27 0.04 0.88 17.67
rain_days 4 -25.78 3.29 0.04 0.92 17.66
minT_totalR 5 -25.04 4.03 0.03 0.94 18.72
minT_maxW 5 -23.83 5.24 0.01 0.96 18.11
maxT_rainDays 5 -23.78 5.29 0.01 0.97 18.09
hail_maxW 5 -23.12 5.95 0.01 0.98 17.76
hail_rainDays 5 -23.03 6.04 0.01 0.99 17.72
lagR_x_maxT 6 -21.36 7.71 0.00 0.99 18.43
minT_x_meanWind 6 -21.30 7.77 0.00 1.00 18.40
minT_x_rainDays 6 -21.04 8.03 0.00 1.00 18.27Check residuals
# identify top-ranked model (no year)
best_no_year_name <- aic_effort_no_year$Modnames[1]
best_no_year_E <- cand_effort_no_year[[best_no_year_name]]
par(mfrow = c(2,2))
plot(best_no_year_E)
# identify top-ranked model (with year)
best_with_year_name <- aic_effort_with_year$Modnames[1]
best_with_year_E <- cand_effort_with_year[[best_with_year_name]]
par(mfrow = c(2,2))
plot(best_with_year_E)
# year-only baseline and best with-year model
summary(m_year)$adj.r.squared[1] 0.758928summary(best_with_year_E)$adj.r.squared[1] 0.7754688# best no-year model
summary(best_no_year_E)$adj.r.squared[1] 0.1105514Figures
ggplot(effort_weather, aes(x = year, y = active_nests)) +
geom_point(size = 2, colour = "grey30") +
geom_smooth(method = "lm", se = TRUE, colour = "black") +
scale_x_continuous(breaks = seq(min(effort_weather$year), max(effort_weather$year), by = 5)) +
labs(
x = "Year",
y = "Number of active nests",
title = "Annual breeding effort at Dronfield Nature Reserve"
) +
theme_bw() +
theme(panel.grid = element_blank())`geom_smooth()` using formula = 'y ~ x'
prep_aic_plot <- function(aic_obj, set_label){
df <- as.data.frame(aic_obj)
# model names
if ("Modnames" %in% names(df)) df$model <- df$Modnames
if (!"model" %in% names(df)) df$model <- rownames(df)
df |>
dplyr::mutate(set = set_label) |>
dplyr::rename(delta = Delta_AICc, weight = AICcWt) |>
dplyr::select(set, model, delta, weight)
}
aic_plot_df <- dplyr::bind_rows(
prep_aic_plot(aic_effort_with_year, "With year"),
prep_aic_plot(aic_effort_no_year, "Without year")
) |>
dplyr::filter(delta <= 3) |>
dplyr::group_by(set) |>
dplyr::mutate(model = forcats::fct_reorder(model, -delta)) |>
dplyr::ungroup()
ggplot(aic_plot_df, aes(x = delta, y = model)) +
geom_vline(xintercept = 2, linetype = 2) +
geom_point(aes(size = weight), colour = "grey20") +
facet_wrap(~set, scales = "free_y") +
scale_size_continuous(range = c(2, 8)) +
labs(
x = expression(Delta*"AICc (lower is better)"),
y = NULL,
title = "Support for candidate models with and without year",
subtitle = "Only models with ΔAICc ≤ 3 are shown; point size indicates AICc weight"
) +
theme_bw() +
theme(panel.grid = element_blank())
# Table 1(weather-only)
aic_df <- as.data.frame(aic_effort_no_year)
# model ids
if ("Modnames" %in% names(aic_df)) {
aic_df <- aic_df |> dplyr::mutate(model_id = Modnames)
} else {
aic_df <- tibble::rownames_to_column(aic_df, "model_id")
}
table1 <- aic_df |>
dplyr::rename(
K = K,
AICc = AICc,
`ΔAICc` = Delta_AICc,
w = AICcWt
) |>
dplyr::select(model_id, K, AICc, `ΔAICc`, w) |>
dplyr::arrange(`ΔAICc`)
# add Adj R2 + slope(SE)
get_stats <- function(mod){
adjr2 <- summary(mod)$adj.r.squared
if (length(coef(mod)) == 2) {
cf <- summary(mod)$coefficients
slope_se <- paste0(round(cf[2,1], 3), " (", round(cf[2,2], 3), ")")
} else {
slope_se <- ""
}
tibble::tibble(`Adj. R²` = adjr2, `slope (SE)` = slope_se)
}
stats_df <- lapply(table1$model_id, function(mn) get_stats(cand_effort_no_year[[mn]])) |>
dplyr::bind_rows() |>
dplyr::mutate(model_id = table1$model_id)
table1 <- table1 |>
dplyr::left_join(stats_df, by = "model_id") |>
dplyr::mutate(
AICc = round(AICc, 2),
`ΔAICc` = round(`ΔAICc`, 2),
w = round(w, 3),
`Adj. R²` = round(`Adj. R²`, 3)
)
# grouping adn neat names
group_order <- c("Baseline","Rain","Temperature","Wind","Hail","Joint effects","Interactions")
model_group <- function(m){
dplyr::case_when(
m == "null" ~ "Baseline",
m %in% c("total_rain","rain_days","max_rain","lag_total_rain") ~ "Rain",
m %in% c("mean_temp","max_temp","min_temp") ~ "Temperature",
m %in% c("mean_wind","max_wind") ~ "Wind",
m %in% c("hail_days") ~ "Hail",
m %in% c("minT_totalR","minT_maxR","minT_maxW","hail_maxW","maxR_maxW","maxT_rainDays","hail_rainDays") ~ "Joint effects",
m %in% c("minT_x_rainDays","minT_x_meanWind","lagR_x_maxT") ~ "Interactions",
TRUE ~ "Joint effects"
)
}
pretty_model <- function(x){
dplyr::recode(
x,
"null" = "Intercept only",
"total_rain" = "Total rainfall",
"rain_days" = "Rainy days",
"max_rain" = "Maximum daily rainfall",
"lag_total_rain" = "Previous year's total rainfall",
"mean_temp" = "Mean temperature",
"max_temp" = "Maximum temperature",
"min_temp" = "Minimum temperature",
"mean_wind" = "Mean wind speed",
"max_wind" = "Maximum wind speed",
"hail_days" = "Hail days",
"minT_totalR" = "Min temperature + total rainfall",
"minT_maxR" = "Min temperature + max daily rainfall",
"minT_maxW" = "Min temperature + max wind speed",
"hail_maxW" = "Hail days + max wind speed",
"maxR_maxW" = "Max daily rainfall + max wind speed",
"maxT_rainDays" = "Max temperature + rainy days",
"hail_rainDays" = "Hail days + rainy days",
"minT_x_rainDays" = "Min temperature × rainy days",
"minT_x_meanWind" = "Min temperature × mean wind speed",
"lagR_x_maxT" = "Prev. rainfall × max temperature",
.default = x
)
}
table1_clean <- table1 |>
dplyr::mutate(
Group = model_group(model_id),
Model = pretty_model(model_id),
Group = factor(Group, levels = group_order)
) |>
dplyr::arrange(Group, `ΔAICc`) |>
dplyr::select(Group, Model, K, AICc, `ΔAICc`, w, `Adj. R²`, `slope (SE)`)
#print
group_sizes <- table1_clean |>
dplyr::count(Group, name = "n") |>
dplyr::mutate(
start = cumsum(dplyr::lag(n, default = 0)) + 1,
end = cumsum(n)
)
tbl <- knitr::kable(
table1_clean |> dplyr::select(-Group),
caption = "Table 1. Weather-only candidate models explaining annual breeding effort (log number of active nests) at Dronfield. K = number of parameters; AICc = small-sample AIC; ΔAICc = difference from the best model; w = Akaike weight. Adjusted R² is shown for all models; slope (SE) is reported only for single-predictor models.",
align = "l"
) |>
kableExtra::kable_styling(full_width = FALSE)
for(i in seq_len(nrow(group_sizes))){
tbl <- tbl |>
kableExtra::pack_rows(as.character(group_sizes$Group[i]),
group_sizes$start[i],
group_sizes$end[i])
}
tbl| Model | K | AICc | ΔAICc | w | Adj. R² | slope (SE) |
|---|---|---|---|---|---|---|
| Baseline | ||||||
| Intercept only | 2 | 14.26 | 2.22 | 0.057 | 0.000 | |
| Rain | ||||||
| Total rainfall | 3 | 12.83 | 0.79 | 0.117 | 0.088 | -0.001 (0) |
| Previous year's total rainfall | 3 | 15.07 | 3.03 | 0.038 | 0.019 | 0 (0) |
| Maximum daily rainfall | 3 | 15.74 | 3.70 | 0.027 | -0.002 | -0.003 (0.003) |
| Rainy days | 3 | 16.72 | 4.68 | 0.017 | -0.034 | 0 (0.003) |
| Temperature | ||||||
| Maximum temperature | 3 | 12.77 | 0.74 | 0.120 | 0.089 | 0.046 (0.023) |
| Mean temperature | 3 | 16.46 | 4.42 | 0.019 | -0.026 | 0.031 (0.062) |
| Minimum temperature | 3 | 16.61 | 4.58 | 0.018 | -0.031 | 0.011 (0.036) |
| Wind | ||||||
| Maximum wind speed | 3 | 14.39 | 2.35 | 0.054 | 0.041 | -0.016 (0.011) |
| Mean wind speed | 3 | 14.84 | 2.81 | 0.043 | 0.026 | -0.18 (0.134) |
| Hail | ||||||
| Hail days | 3 | 12.04 | 0.00 | 0.174 | 0.111 | 0.076 (0.035) |
| Joint effects | ||||||
| Hail days + max wind speed | 4 | 13.16 | 1.12 | 0.099 | 0.123 | |
| Max temperature + rainy days | 4 | 14.37 | 2.33 | 0.054 | 0.088 | |
| Hail days + rainy days | 4 | 14.44 | 2.40 | 0.052 | 0.086 | |
| Min temperature + total rainfall | 4 | 15.45 | 3.41 | 0.032 | 0.056 | |
| Max daily rainfall + max wind speed | 4 | 15.84 | 3.81 | 0.026 | 0.044 | |
| Min temperature + max wind speed | 4 | 17.03 | 5.00 | 0.014 | 0.006 | |
| Min temperature + max daily rainfall | 4 | 18.24 | 6.20 | 0.008 | -0.033 | |
| Interactions | ||||||
| Prev. rainfall × max temperature | 5 | 15.71 | 3.67 | 0.028 | 0.100 | |
| Min temperature × mean wind speed | 5 | 20.27 | 8.23 | 0.003 | -0.043 | |
| Min temperature × rainy days | 5 | 21.95 | 9.92 | 0.001 | -0.101 | |
digits_aic <- 2
digits_r2 <- 3
digits_slope <- 3
group_order <- c("Baseline", "Rain", "Temperature", "Wind", "Hail", "Joint effects", "Interactions")
model_group_with_year <- function(m){
dplyr::case_when(
m %in% c("year") ~ "Baseline",
m %in% c("total_rain","rain_days","max_rain","lag_total_rain") ~ "Rain",
m %in% c("mean_temp","max_temp","min_temp") ~ "Temperature",
m %in% c("mean_wind","max_wind") ~ "Wind",
m %in% c("hail_days") ~ "Hail",
m %in% c("minT_totalR","minT_maxR","minT_maxW","hail_maxW","maxR_maxW","maxT_rainDays","hail_rainDays") ~ "Joint effects",
m %in% c("minT_x_rainDays","minT_x_meanWind","lagR_x_maxT") ~ "Interactions",
TRUE ~ "Joint effects"
)
}
pretty_model_with_year <- function(x){
dplyr::recode(
x,
"year" = "Year (baseline)",
"total_rain" = "Total rainfall",
"rain_days" = "Rainy days",
"max_rain" = "Maximum daily rainfall",
"lag_total_rain" = "Previous year's total rainfall",
"mean_temp" = "Mean temperature",
"max_temp" = "Maximum temperature",
"min_temp" = "Minimum temperature",
"mean_wind" = "Mean wind speed",
"max_wind" = "Maximum wind speed",
"hail_days" = "Hail days",
"minT_totalR" = "Min temperature + total rainfall",
"minT_maxR" = "Min temperature + max daily rainfall",
"minT_maxW" = "Min temperature + max wind speed",
"hail_maxW" = "Hail days + max wind speed",
"maxR_maxW" = "Max daily rainfall + max wind speed",
"maxT_rainDays" = "Max temperature + rainy days",
"hail_rainDays" = "Hail days + rainy days",
"minT_x_rainDays" = "Min temperature × rainy days",
"minT_x_meanWind" = "Min temperature × mean wind speed",
"lagR_x_maxT" = "Prev. rainfall × max temperature",
.default = x
)
}
aic_df2 <- as.data.frame(aic_effort_with_year)
if ("Modnames" %in% names(aic_df2)) {
aic_df2 <- aic_df2 |> dplyr::mutate(model = Modnames)
} else {
aic_df2 <- tibble::rownames_to_column(aic_df2, "model")
}
aic_df2 <- aic_df2 |>
dplyr::rename(
K = K,
AICc = AICc,
dAICc = Delta_AICc,
w = AICcWt
) |>
dplyr::select(model, K, AICc, dAICc, w) |>
dplyr::arrange(dAICc)
get_stats_with_year <- function(mod){
adjr2 <- summary(mod)$adj.r.squared
if (length(coef(mod)) == 3) { # intercept + year + ONE weather predictor
cf <- summary(mod)$coefficients
slope <- unname(cf[3, 1])
se <- unname(cf[3, 2])
slope_se <- paste0(
round(slope, digits_slope),
" (", round(se, digits_slope), ")"
)
} else {
slope_se <- ""
}
tibble::tibble(adj_R2 = adjr2, slope_SE = slope_se)
}
stats_df2 <- lapply(aic_df2$model, function(mn) get_stats_with_year(cand_effort_with_year[[mn]])) |>
dplyr::bind_rows() |>
dplyr::mutate(model = aic_df2$model)
table2 <- aic_df2 |>
dplyr::left_join(stats_df2, by = "model") |>
dplyr::mutate(
AICc = round(AICc, digits_aic),
dAICc = round(dAICc, digits_aic),
w = round(w, 3),
adj_R2 = round(adj_R2, digits_r2)
) |>
dplyr::mutate(
model_id = model,
Group = model_group_with_year(model_id),
Model = pretty_model_with_year(model_id),
Group = factor(Group, levels = group_order)
) |>
dplyr::arrange(Group, dAICc) |>
dplyr::rename(
`ΔAICc` = dAICc,
`Adj. R²` = adj_R2,
`slope (SE)` = slope_SE
) |>
dplyr::select(Group, Model, K, AICc, `ΔAICc`, w, `Adj. R²`, `slope (SE)`)
group_sizes2 <- table2 |>
dplyr::count(Group, name = "n") |>
dplyr::filter(n > 0) |>
dplyr::mutate(
start = cumsum(dplyr::lag(n, default = 0)) + 1,
end = cumsum(n)
)
tbl2 <- knitr::kable(
table2 |> dplyr::select(-Group),
caption = "Table 2. Candidate models including year explaining annual breeding effort (log number of active nests) at Dronfield. K = number of parameters; AICc = small-sample AIC; ΔAICc = difference from the best model; w = Akaike weight. Adjusted R² is shown for all models; slope (SE) is reported for each model (non-intercept terms shown as estimate (SE)).",
align = "l"
) |>
kableExtra::kable_styling(full_width = FALSE)
for(i in seq_len(nrow(group_sizes2))){
tbl2 <- tbl2 |>
kableExtra::pack_rows(
as.character(group_sizes2$Group[i]),
group_sizes2$start[i],
group_sizes2$end[i]
)
}
tbl2| Model | K | AICc | ΔAICc | w | Adj. R² | slope (SE) |
|---|---|---|---|---|---|---|
| Baseline | ||||||
| Year (baseline) | 3 | -28.43 | 0.64 | 0.140 | 0.759 | |
| Rain | ||||||
| Maximum daily rainfall | 4 | -29.07 | 0.00 | 0.193 | 0.775 | -0.002 (0.001) |
| Total rainfall | 4 | -27.49 | 1.58 | 0.088 | 0.764 | 0 (0) |
| Previous year's total rainfall | 4 | -25.81 | 3.27 | 0.038 | 0.751 | 0 (0) |
| Rainy days | 4 | -25.78 | 3.29 | 0.037 | 0.750 | 0 (0.002) |
| Temperature | ||||||
| Minimum temperature | 4 | -26.68 | 2.39 | 0.058 | 0.757 | 0.016 (0.017) |
| Mean temperature | 4 | -26.59 | 2.48 | 0.056 | 0.757 | 0.026 (0.03) |
| Maximum temperature | 4 | -26.42 | 2.65 | 0.051 | 0.755 | 0.01 (0.013) |
| Wind | ||||||
| Maximum wind speed | 4 | -25.90 | 3.18 | 0.039 | 0.751 | -0.002 (0.006) |
| Mean wind speed | 4 | -25.87 | 3.21 | 0.039 | 0.751 | -0.019 (0.07) |
| Hail | ||||||
| Hail days | 4 | -25.88 | 3.19 | 0.039 | 0.751 | 0.006 (0.02) |
| Joint effects | ||||||
| Min temperature + max daily rainfall | 5 | -27.44 | 1.63 | 0.085 | 0.776 | |
| Max daily rainfall + max wind speed | 5 | -26.41 | 2.66 | 0.051 | 0.769 | |
| Min temperature + total rainfall | 5 | -25.04 | 4.03 | 0.026 | 0.758 | |
| Min temperature + max wind speed | 5 | -23.83 | 5.24 | 0.014 | 0.749 | |
| Max temperature + rainy days | 5 | -23.78 | 5.29 | 0.014 | 0.748 | |
| Hail days + max wind speed | 5 | -23.12 | 5.95 | 0.010 | 0.743 | |
| Hail days + rainy days | 5 | -23.03 | 6.04 | 0.009 | 0.742 | |
| Interactions | ||||||
| Prev. rainfall × max temperature | 6 | -21.36 | 7.71 | 0.004 | 0.744 | |
| Min temperature × mean wind speed | 6 | -21.30 | 7.77 | 0.004 | 0.744 | |
| Min temperature × rainy days | 6 | -21.04 | 8.03 | 0.003 | 0.741 | |
Objective 2: Breeding success
objective2_table <- tribble(
~`Hypothesis / justification`, ~Model,
# ---- Temperature ----
"**Temp**", "",
"Extreme heat events cause chick heat stress and dehydration, reducing fledging success",
"success ~ max_temp + (1 | year)",
"Cold extremes increase thermoregulatory demand in chicks, increasing mortality risk",
"success ~ min_temp + (1 | year)",
# ---- Rain ----
"**Rain**", "",
"Prolonged rainfall reduces adult foraging efficiency, limiting food delivery to chicks and reducing breeding success",
"success ~ rain_days + (1 | year)",
"Intense rainfall events cause acute nest disturbance and chick exposure, reducing breeding success",
"success ~ max_rain + (1 | year)",
"Overall wet years reduce provisioning efficiency across the breeding season, lowering breeding success",
"success ~ total_rain + (1 | year)",
# ---- Wind ----
"**Wind**", "",
"Extreme winds increase nest exposure and disrupt provisioning, reducing chick survival",
"success ~ max_wind + (1 | year)",
# ---- Hail ----
"**Hail**", "",
"Severe hail events increase nest destruction and chick mortality, reducing breeding success",
"success ~ hail_days + (1 | year)",
# ---- Joint Effects ----
"**Joint effects**", "",
"Severe storm exposure increases nest destruction and chick mortality, reducing breeding success",
"success ~ hail_days + max_wind + (1 | year)",
"Prolonged wet conditions combined with storm events increase chick exposure and limit provisioning, reducing breeding success",
"success ~ rain_days + hail_days + (1 | year)",
"Heat stress combined with prolonged wet conditions reduces chick thermoregulation and food delivery, lowering breeding success",
"success ~ max_temp + rain_days + (1 | year)",
# ---- Interactive Effects ----
"**Interactive effects**", "",
"The negative effect of heat stress on breeding success is amplified during persistently wet conditions that limit adult foraging",
"success ~ max_temp * rain_days + (1 | year)",
"Cold stress on chicks is amplified during extreme wind events due to increased exposure and heat loss, reducing breeding success",
"success ~ min_temp * max_wind + (1 | year)"
)
kable(
objective2_table,
align = c("l", "l"),
caption = "Objective 2 candidate model set for breeding success. All models are fitted as binomial GLMMs: cbind(success_nests, failed_nests) with a random intercept for year (1 | year)"
)| Hypothesis / justification | Model |
|---|---|
| Temp | |
| Extreme heat events cause chick heat stress and dehydration, reducing fledging success | success ~ max_temp + (1 | year) |
| Cold extremes increase thermoregulatory demand in chicks, increasing mortality risk | success ~ min_temp + (1 | year) |
| Rain | |
| Prolonged rainfall reduces adult foraging efficiency, limiting food delivery to chicks and reducing breeding success | success ~ rain_days + (1 | year) |
| Intense rainfall events cause acute nest disturbance and chick exposure, reducing breeding success | success ~ max_rain + (1 | year) |
| Overall wet years reduce provisioning efficiency across the breeding season, lowering breeding success | success ~ total_rain + (1 | year) |
| Wind | |
| Extreme winds increase nest exposure and disrupt provisioning, reducing chick survival | success ~ max_wind + (1 | year) |
| Hail | |
| Severe hail events increase nest destruction and chick mortality, reducing breeding success | success ~ hail_days + (1 | year) |
| Joint effects | |
| Severe storm exposure increases nest destruction and chick mortality, reducing breeding success | success ~ hail_days + max_wind + (1 | year) |
| Prolonged wet conditions combined with storm events increase chick exposure and limit provisioning, reducing breeding success | success ~ rain_days + hail_days + (1 | year) |
| Heat stress combined with prolonged wet conditions reduces chick thermoregulation and food delivery, lowering breeding success | success ~ max_temp + rain_days + (1 | year) |
| Interactive effects | |
| The negative effect of heat stress on breeding success is amplified during persistently wet conditions that limit adult foraging | success ~ max_temp * rain_days + (1 | year) |
| Cold stress on chicks is amplified during extreme wind events due to increased exposure and heat loss, reducing breeding success | success ~ min_temp * max_wind + (1 | year) |
success ~ 1 + (1 | year)
library(lme4)
success_df <- success_df |> mutate(year = factor(year))
m0 <- glmer(
cbind(success_nests, failed_nests) ~ 1 + (1 | year),
family = binomial,
data = success_df
)
summary(m0)Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: cbind(success_nests, failed_nests) ~ 1 + (1 | year)
Data: success_df
AIC BIC logLik -2*log(L) df.resid
229.0 232.0 -112.5 225.0 30
Scaled residuals:
Min 1Q Median 3Q Max
-1.25016 -0.28770 0.00005 0.48287 1.24098
Random effects:
Groups Name Variance Std.Dev.
year (Intercept) 0.1292 0.3594
Number of obs: 32, groups: year, 32
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.21944 0.07595 2.889 0.00386 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(lme4)
success_df <- success_df |> dplyr::mutate(year = factor(year))
# Baseline (null) model
m0_success <- glmer(
cbind(success_nests, failed_nests) ~ 1 + (1 | year),
family = binomial,
data = success_df
)
# Single predictor models
m_succ_maxT <- glmer(
cbind(success_nests, failed_nests) ~ z_max_temp + (1 | year),
family = binomial, data = success_df
)
m_succ_minT <- glmer(
cbind(success_nests, failed_nests) ~ z_min_temp + (1 | year),
family = binomial, data = success_df
)
m_succ_rainD <- glmer(
cbind(success_nests, failed_nests) ~ z_rain_days + (1 | year),
family = binomial, data = success_df
)
m_succ_maxR <- glmer(
cbind(success_nests, failed_nests) ~ z_max_rain + (1 | year),
family = binomial, data = success_df
)
m_succ_totR <- glmer(
cbind(success_nests, failed_nests) ~ z_total_rain + (1 | year),
family = binomial, data = success_df
)
m_succ_maxW <- glmer(
cbind(success_nests, failed_nests) ~ z_max_wind + (1 | year),
family = binomial, data = success_df
)
m_succ_hail <- glmer(
cbind(success_nests, failed_nests) ~ z_hail_days + (1 | year),
family = binomial, data = success_df
)
# Joint effects
m_succ_hailW <- glmer(
cbind(success_nests, failed_nests) ~ z_hail_days + z_max_wind + (1 | year),
family = binomial, data = success_df
)
m_succ_rainH <- glmer(
cbind(success_nests, failed_nests) ~ z_rain_days + z_hail_days + (1 | year),
family = binomial, data = success_df
)
m_succ_heatR <- glmer(
cbind(success_nests, failed_nests) ~ z_max_temp + z_rain_days + (1 | year),
family = binomial, data = success_df
)
# Interactions
m_succ_int_heat_rain <- glmer(
cbind(success_nests, failed_nests) ~ z_max_temp * z_rain_days + (1 | year),
family = binomial, data = success_df
)
m_succ_int_cold_wind <- glmer(
cbind(success_nests, failed_nests) ~ z_min_temp * z_max_wind + (1 | year),
family = binomial, data = success_df
)# Candidate set + AIC
library(AICcmodavg)
cand_success <- list(
null_re = m0_success,
max_temp = m_succ_maxT,
min_temp = m_succ_minT,
rain_days = m_succ_rainD,
max_rain = m_succ_maxR,
total_rain = m_succ_totR,
max_wind = m_succ_maxW,
hail_days = m_succ_hail,
hail_maxW = m_succ_hailW,
rain_hail = m_succ_rainH,
maxT_rainD = m_succ_heatR,
int_heatRain = m_succ_int_heat_rain,
int_coldWind = m_succ_int_cold_wind
)
aic_success <- aictab(
cand.set = cand_success,
modnames = names(cand_success)
)
aic_success
Model selection based on AICc:
K AICc Delta_AICc AICcWt Cum.Wt LL
min_temp 3 226.93 0.00 0.43 0.43 -110.03
null_re 2 229.45 2.52 0.12 0.56 -112.52
hail_days 3 230.45 3.53 0.07 0.63 -111.80
max_temp 3 230.85 3.92 0.06 0.69 -111.99
max_rain 3 231.01 4.09 0.06 0.75 -112.08
rain_days 3 231.61 4.68 0.04 0.79 -112.38
total_rain 3 231.65 4.72 0.04 0.83 -112.40
max_wind 3 231.73 4.81 0.04 0.87 -112.44
int_coldWind 5 232.04 5.11 0.03 0.91 -109.87
maxT_rainD 4 232.10 5.18 0.03 0.94 -111.31
rain_hail 4 232.30 5.38 0.03 0.97 -111.41
hail_maxW 4 232.67 5.74 0.02 0.99 -111.59
int_heatRain 5 234.93 8.01 0.01 1.00 -111.31vc_year <- function(mod) {
vc <- as.data.frame(VarCorr(mod))
vc$vcov[vc$grp == "year"][1]
}
v0 <- vc_year(m0_success)
v_minT <- vc_year(m_succ_minT)
prop_between_year_explained <- (v0 - v_minT) / v0
prop_between_year_explained[1] 0.2098479Figures and tables
library(dplyr)
library(tibble)
library(purrr)
library(knitr)
library(kableExtra)
Attaching package: 'kableExtra'The following object is masked from 'package:dplyr':
group_rowslibrary(lme4)
aic_df <- as.data.frame(aic_success)
if ("Modnames" %in% names(aic_df)) {
aic_df <- aic_df |> mutate(model_id = Modnames)
} else {
aic_df <- tibble::rownames_to_column(aic_df, "model_id")
}
table_success <- aic_df |>
rename(
K = K,
AICc = AICc,
`ΔAICc` = Delta_AICc,
w = AICcWt
) |>
dplyr::select(model_id, K, AICc, `ΔAICc`, w) |>
arrange(`ΔAICc`)
vc_year <- function(mod) {
vc <- as.data.frame(VarCorr(mod))
vc$vcov[vc$grp == "year"][1]
}
v0 <- vc_year(m0_success)
var_expl_tbl <- tibble(
model_id = names(cand_success),
`Variance explained` = map_dbl(cand_success, ~ {
v <- vc_year(.x)
(v0 - v) / v0
})
) |>
mutate(`Variance explained` = round(`Variance explained`, 3))
table_success <- table_success |>
left_join(var_expl_tbl, by = "model_id")
get_slope_se <- function(mod) {
cf <- summary(mod)$coefficients
if (nrow(cf) == 2) {
paste0(round(cf[2, 1], 3), " (", round(cf[2, 2], 3), ")")
} else {
""
}
}
slope_tbl <- tibble(
model_id = table_success$model_id,
`slope (SE)` = map_chr(table_success$model_id, ~ get_slope_se(cand_success[[.x]]))
)
table_success <- table_success |>
left_join(slope_tbl, by = "model_id") |>
mutate(
AICc = round(AICc, 2),
`ΔAICc` = round(`ΔAICc`, 2),
w = round(w, 3)
)
group_order <- c("Baseline", "Temperature", "Rain", "Wind", "Hail", "Joint effects", "Interactions")
model_group_succ <- function(m) {
case_when(
m %in% c("null_re", "null", "null_success") ~ "Baseline",
m %in% c("min_temp", "max_temp") ~ "Temperature",
m %in% c("rain_days", "max_rain", "total_rain") ~ "Rain",
m %in% c("max_wind") ~ "Wind",
m %in% c("hail_days") ~ "Hail",
m %in% c("hail_maxW", "rain_hail", "maxT_rainD") ~ "Joint effects",
m %in% c("int_heatRain", "int_coldWind") ~ "Interactions",
TRUE ~ "Joint effects"
)
}
neat_succ <- function(x) {
recode(
x,
"null_re" = "Intercept only",
"min_temp" = "Minimum temperature",
"max_temp" = "Maximum temperature",
"rain_days" = "Rainy days",
"max_rain" = "Maximum daily rainfall",
"total_rain" = "Total rainfall",
"max_wind" = "Maximum wind speed",
"hail_days" = "Hail days",
"hail_maxW" = "Hail days + max wind speed",
"rain_hail" = "Rainy days + hail days",
"maxT_rainD" = "Max temperature + rainy days",
"int_heatRain" = "Max temperature × rainy days",
"int_coldWind" = "Min temperature × max wind speed",
.default = x
)
}
table_success_clean <- table_success |>
mutate(
Group = model_group_succ(model_id),
Model = neat_succ(model_id),
Group = factor(Group, levels = group_order)
) |>
arrange(Group, `ΔAICc`) |>
dplyr::select(Group, Model, K, AICc, `ΔAICc`, w, `Variance explained`, `slope (SE)`)
group_sizes <- table_success_clean |>
count(Group, name = "n") |>
filter(n > 0) |>
mutate(
start = cumsum(lag(n, default = 0)) + 1,
end = cumsum(n)
)
tbl_succ <- kable(
table_success_clean |> dplyr::select(-Group),
caption = "Table 3:Model selection results for breeding success (binomial GLMMs with a random intercept for year). Variance explained is the proportional reduction in the estimated between-year variance relative to the null model. Slope (SE) is shown only for single-predictor models.",
align = "l"
) |>
kable_styling(full_width = FALSE)
for (i in seq_len(nrow(group_sizes))) {
tbl_succ <- tbl_succ |>
pack_rows(
as.character(group_sizes$Group[i]),
group_sizes$start[i],
group_sizes$end[i]
)
}
tbl_succ| Model | K | AICc | ΔAICc | w | Variance explained | slope (SE) |
|---|---|---|---|---|---|---|
| Baseline | ||||||
| Intercept only | 2 | 229.45 | 2.52 | 0.123 | 0.000 | |
| Temperature | ||||||
| Minimum temperature | 3 | 226.93 | 0.00 | 0.435 | 0.210 | 0.167 (0.072) |
| Maximum temperature | 3 | 230.85 | 3.92 | 0.061 | 0.042 | -0.077 (0.075) |
| Rain | ||||||
| Maximum daily rainfall | 3 | 231.01 | 4.09 | 0.056 | 0.034 | 0.075 (0.079) |
| Rainy days | 3 | 231.61 | 4.68 | 0.042 | 0.014 | -0.041 (0.076) |
| Total rainfall | 3 | 231.65 | 4.72 | 0.041 | 0.008 | 0.038 (0.076) |
| Wind | ||||||
| Maximum wind speed | 3 | 231.73 | 4.81 | 0.039 | 0.003 | -0.031 (0.077) |
| Hail | ||||||
| Hail days | 3 | 230.45 | 3.53 | 0.074 | 0.052 | -0.088 (0.073) |
| Joint effects | ||||||
| Max temperature + rainy days | 4 | 232.10 | 5.18 | 0.033 | 0.104 | |
| Rainy days + hail days | 4 | 232.30 | 5.38 | 0.030 | 0.087 | |
| Hail days + max wind speed | 4 | 232.67 | 5.74 | 0.025 | 0.062 | |
| Interactions | ||||||
| Min temperature × max wind speed | 5 | 232.04 | 5.11 | 0.034 | 0.225 | |
| Max temperature × rainy days | 5 | 234.93 | 8.01 | 0.008 | 0.104 | |
minT_model <- m_succ_minT
success_df_plot <- success_df |>
mutate(
success_rate = success_nests / (success_nests + failed_nests),
year_num = as.numeric(as.character(year))
)
yrs <- sort(unique(success_df_plot$year_num[is.finite(success_df_plot$year_num)]))
breaks_5y <- yrs[seq(1, length(yrs), by = 5)] #
ggplot(success_df_plot, aes(x = year_num, y = success_rate)) +
geom_point(size = 2, colour = "grey30") +
geom_smooth(
method = "lm",
se = TRUE,
colour = "grey30",
fill = "grey85"
) +
scale_y_continuous(limits = c(0, 1)) +
scale_x_continuous(breaks = breaks_5y) +
labs(
x = "Year",
y = "Breeding success",
title = "Annual breeding success (proportion of successful nests) of White-backed Vultures"
) +
theme_bw() +
theme(panel.grid = element_blank())`geom_smooth()` using formula = 'y ~ x'
obs_df <- success_df |>
mutate(
success_rate = success_nests / (success_nests + failed_nests),
se = sqrt(success_rate * (1 - success_rate) /
(success_nests + failed_nests))
)minT_model <- m_succ_minT
pred_df <- success_df |>
mutate(
pred = predict(
minT_model,
newdata = success_df,
type = "response",
re.form = NA
)
)obs_df <- success_df |>
mutate(
year_num = as.numeric(gsub("[^0-9]", "", as.character(year))),
success_rate = success_nests / (success_nests + failed_nests),
se = sqrt(success_rate * (1 - success_rate) / (success_nests + failed_nests))
)
minT_model <- m_succ_minT
pred_df <- success_df |>
mutate(
year_num = as.numeric(gsub("[^0-9]", "", as.character(year))),
pred = predict(minT_model, newdata = success_df, type = "response", re.form = NA)
)
ggplot() +
geom_point(
data = obs_df,
aes(x = year_num, y = success_rate),
size = 2, shape = 21, fill = "white", colour = "grey30"
) +
geom_errorbar(
data = obs_df,
aes(x = year_num,
ymin = success_rate - se,
ymax = success_rate + se),
width = 0, colour = "grey30"
) +
geom_line(
data = pred_df,
aes(x = year_num, y = pred),
linewidth = 1, colour = "black"
) +
scale_y_continuous(limits = c(0, 1)) +
scale_x_continuous(breaks = seq(1995, 2025, by = 5)) +
labs(
x = "Year",
y = "Breeding success",
title = "Observed and temperature-predicted breeding success",
subtitle = "Points show annual observed success (± SE); line shows predictions from the minimum temperature model"
) +
theme_bw() +
theme(panel.grid = element_blank())
figure shows that years with warmer minimum temperatures tend to have slightly higher breeding success, but temperature alone does not explain most of the year-to-year variation. Given the temperature in each year, this is what the model predicts breeding success should be
Objective 3
hist(timing_weather_60$lay_DOY)
qqnorm(timing_weather_60$lay_DOY)
qqline(timing_weather_60$lay_DOY)
objective3_table <- tribble(
~`Hypothesis / justification`, ~Model,
# ---- Temperature (pre-laying cues) ----
"**Temp (pre-laying)**", "",
"Short-term warm conditions prior to laying may act as a cue for the initiation of breeding.",
"lay_DOY ~ mean_temp_60",
"Cold pre-laying conditions delay breeding due to increased energetic costs",
"lay_DOY ~ min_temp_60",
# ---- Rain (pre-laying cues) ----
"**Rain (pre-laying)**", "",
"Prolonged wet conditions prior to laying delay breeding by limiting foraging efficiency",
"lay_DOY ~ rain_days_60",
"Higher rainfall prior to laying signals improved food availability and advances breeding timing",
"lay_DOY ~ total_rain_60",
# ---- Wind (pre-laying flight conditions) ----
"**Wind (pre-laying)**", "",
"Persistent windy conditions prior to laying increase flight costs and delay breeding timing",
"lay_DOY ~ mean_wind_60",
# ---- Lagged effects (carry-over cues) ----
"**Lagged effects (previous year)**", "",
"Higher rainfall in the previous year improves food availability and adult condition entering the breeding season, allowing vultures to initiate breeding earlier",
"lay_DOY ~ lag_total_rain",
"Prolonged wet conditions in the previous year may influence adult condition and carry over to affect laying timing",
"lay_DOY ~ lag_rain_days",
"Thermal conditions in the previous year influence physiological state entering the breeding season",
"lay_DOY ~ lag_mean_temp",
# ---- Joint effects ----
"**Joint effects**", "",
"Temperature and rainfall cues prior to laying jointly influence breeding timing",
"lay_DOY ~ mean_temp_60 + rain_days_60",
"Carry-over effects from the previous year interact with current pre-laying rainfall cues",
"lay_DOY ~ lag_total_rain + rain_days_60",
"Previous-year rainfall and current temperature jointly influence breeding timing",
"lay_DOY ~ lag_total_rain + mean_temp_60",
# ---- Interactive effects ----
"**Interactive effects**", "",
"The effect of temperature cues on breeding timing depends on pre-laying rainfall conditions",
"lay_DOY ~ mean_temp_60 * rain_days_60",
"Carry-over effects from the previous year modify responses to current pre-laying temperature cues",
"lay_DOY ~ lag_total_rain * mean_temp_60"
)
kable(
objective3_table,
align = c("l", "l"),
caption = "Objective 3 candidate model set for breeding timing (laying date, day of year)."
)| Hypothesis / justification | Model |
|---|---|
| Temp (pre-laying) | |
| Short-term warm conditions prior to laying may act as a cue for the initiation of breeding. | lay_DOY ~ mean_temp_60 |
| Cold pre-laying conditions delay breeding due to increased energetic costs | lay_DOY ~ min_temp_60 |
| Rain (pre-laying) | |
| Prolonged wet conditions prior to laying delay breeding by limiting foraging efficiency | lay_DOY ~ rain_days_60 |
| Higher rainfall prior to laying signals improved food availability and advances breeding timing | lay_DOY ~ total_rain_60 |
| Wind (pre-laying) | |
| Persistent windy conditions prior to laying increase flight costs and delay breeding timing | lay_DOY ~ mean_wind_60 |
| Lagged effects (previous year) | |
| Higher rainfall in the previous year improves food availability and adult condition entering the breeding season, allowing vultures to initiate breeding earlier | lay_DOY ~ lag_total_rain |
| Prolonged wet conditions in the previous year may influence adult condition and carry over to affect laying timing | lay_DOY ~ lag_rain_days |
| Thermal conditions in the previous year influence physiological state entering the breeding season | lay_DOY ~ lag_mean_temp |
| Joint effects | |
| Temperature and rainfall cues prior to laying jointly influence breeding timing | lay_DOY ~ mean_temp_60 + rain_days_60 |
| Carry-over effects from the previous year interact with current pre-laying rainfall cues | lay_DOY ~ lag_total_rain + rain_days_60 |
| Previous-year rainfall and current temperature jointly influence breeding timing | lay_DOY ~ lag_total_rain + mean_temp_60 |
| Interactive effects | |
| The effect of temperature cues on breeding timing depends on pre-laying rainfall conditions | lay_DOY ~ mean_temp_60 * rain_days_60 |
| Carry-over effects from the previous year modify responses to current pre-laying temperature cues | lay_DOY ~ lag_total_rain * mean_temp_60 |
# 60 days prior
timing_weather_60 <- timing_weather_60 |>
left_join(
weather_yearly_lag |>
dplyr::select(
YEAR,
lag_total_rain,
lag_rain_days,
lag_mean_temp
),
by = c("year" = "YEAR")
)
# 30 days prior
timing_weather_30 <- timing_weather_30 |>
left_join(
weather_yearly_lag |>
dplyr::select(
YEAR,
lag_total_rain,
lag_rain_days,
lag_mean_temp
),
by = c("year" = "YEAR")
)cand_timing_30 <- list(
meanT_30 = lm(lay_DOY ~ mean_temp_30, data = timing_weather_30),
maxT_30 = lm(lay_DOY ~ max_temp_30, data = timing_weather_30),
minT_30 = lm(lay_DOY ~ min_temp_30, data = timing_weather_30),
rainDays_30 = lm(lay_DOY ~ rain_days_30, data = timing_weather_30),
totalR_30 = lm(lay_DOY ~ total_rain_30, data = timing_weather_30),
meanW_30 = lm(lay_DOY ~ mean_wind_30, data = timing_weather_30),
maxW_30 = lm(lay_DOY ~ max_wind_30, data = timing_weather_30),
# lagged (previous year)
lagTotR = lm(lay_DOY ~ lag_total_rain, data = timing_weather_30),
lagRainD = lm(lay_DOY ~ lag_rain_days, data = timing_weather_30),
lagMeanT = lm(lay_DOY ~ lag_mean_temp, data = timing_weather_30),
# additive
temp_rain = lm(lay_DOY ~ mean_temp_30 + rain_days_30, data = timing_weather_30),
temp_wind = lm(lay_DOY ~ mean_temp_30 + mean_wind_30, data = timing_weather_30),
rain_wind = lm(lay_DOY ~ rain_days_30 + mean_wind_30, data = timing_weather_30),
# interactions
temp_x_rain = lm(lay_DOY ~ mean_temp_30 * rain_days_30, data = timing_weather_30),
minT_x_wind = lm(lay_DOY ~ min_temp_30 * mean_wind_30, data = timing_weather_30),
lagR_x_temp = lm(lay_DOY ~ lag_total_rain * mean_temp_30, data = timing_weather_30)
)aic_timing_30 <- aictab(
cand.set = cand_timing_30,
modnames = names(cand_timing_30)
)cand_timing_60 <- list(
meanT_60 = lm(lay_DOY ~ mean_temp_60, data = timing_weather_60),
maxT_60 = lm(lay_DOY ~ max_temp_60, data = timing_weather_60),
minT_60 = lm(lay_DOY ~ min_temp_60, data = timing_weather_60),
rainDays_60 = lm(lay_DOY ~ rain_days_60, data = timing_weather_60),
totalR_60 = lm(lay_DOY ~ total_rain_60, data = timing_weather_60),
meanW_60 = lm(lay_DOY ~ mean_wind_60, data = timing_weather_60),
maxW_60 = lm(lay_DOY ~ max_wind_60, data = timing_weather_60),
# lagged (previous year)
lagTotR = lm(lay_DOY ~ lag_total_rain, data = timing_weather_60),
lagRainD = lm(lay_DOY ~ lag_rain_days, data = timing_weather_60),
lagMeanT = lm(lay_DOY ~ lag_mean_temp, data = timing_weather_60),
# additive
temp_rain = lm(lay_DOY ~ mean_temp_60 + rain_days_60, data = timing_weather_60),
temp_wind = lm(lay_DOY ~ mean_temp_60 + mean_wind_60, data = timing_weather_60),
rain_wind = lm(lay_DOY ~ rain_days_60 + mean_wind_60, data = timing_weather_60),
# interactions
temp_x_rain = lm(lay_DOY ~ mean_temp_60 * rain_days_60, data = timing_weather_60),
minT_x_wind = lm(lay_DOY ~ min_temp_60 * mean_wind_60, data = timing_weather_60),
lagR_x_temp = lm(lay_DOY ~ lag_total_rain * mean_temp_60, data = timing_weather_60)
)aic_timing_60 <- aictab(
cand.set = cand_timing_60,
modnames = names(cand_timing_60)
)aic_timing_30
Model selection based on AICc:
K AICc Delta_AICc AICcWt Cum.Wt LL
temp_x_rain 5 11336.00 0.00 1 1 -5662.98
temp_rain 4 11349.65 13.65 0 1 -5670.81
lagR_x_temp 5 11373.62 37.62 0 1 -5681.79
minT_x_wind 5 11438.51 102.51 0 1 -5714.24
minT_30 3 11529.94 193.94 0 1 -5761.96
temp_wind 4 11571.91 235.91 0 1 -5781.94
meanT_30 3 11600.33 264.33 0 1 -5797.16
maxT_30 3 11820.41 484.41 0 1 -5907.20
rain_wind 4 12224.48 888.48 0 1 -6108.23
lagTotR 3 12278.29 942.29 0 1 -6136.14
lagRainD 3 12289.62 953.61 0 1 -6141.80
lagMeanT 3 12292.70 956.70 0 1 -6143.34
rainDays_30 3 12360.96 1024.96 0 1 -6177.47
totalR_30 3 12400.54 1064.54 0 1 -6197.26
maxW_30 3 12506.89 1170.89 0 1 -6250.44
meanW_30 3 12516.32 1180.31 0 1 -6255.15summary(cand_timing_30$temp_x_rain)$r.squared[1] 0.5628433aic_timing_60
Model selection based on AICc:
K AICc Delta_AICc AICcWt Cum.Wt LL
temp_x_rain 5 10760.71 0.00 0.9 0.9 -5375.33
temp_rain 4 10765.17 4.46 0.1 1.0 -5378.57
lagR_x_temp 5 10889.01 128.31 0.0 1.0 -5439.49
minT_x_wind 5 10945.22 184.51 0.0 1.0 -5467.59
minT_60 3 10990.38 229.67 0.0 1.0 -5492.18
temp_wind 4 11100.18 339.48 0.0 1.0 -5546.08
meanT_60 3 11112.60 351.89 0.0 1.0 -5553.29
maxT_60 3 11394.57 633.86 0.0 1.0 -5694.28
rain_wind 4 12085.12 1324.41 0.0 1.0 -6038.55
rainDays_60 3 12211.31 1450.61 0.0 1.0 -6102.65
lagTotR 3 12278.29 1517.58 0.0 1.0 -6136.14
totalR_60 3 12282.33 1521.62 0.0 1.0 -6138.15
lagRainD 3 12289.62 1528.91 0.0 1.0 -6141.80
lagMeanT 3 12292.70 1531.99 0.0 1.0 -6143.34
maxW_60 3 12533.72 1773.01 0.0 1.0 -6263.85
meanW_60 3 12538.92 1778.21 0.0 1.0 -6266.45summary(cand_timing_60$temp_x_rain)$r.squared[1] 0.7052922m_year_trend <- lm(lay_DOY ~ year, data = timing_weather_30)
summary(m_year_trend)
Call:
lm(formula = lay_DOY ~ year, data = timing_weather_30)
Residuals:
Min 1Q Median 3Q Max
-145.506 -10.576 -2.616 7.142 142.283
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 51.56562 99.65756 0.517 0.605
year 0.05049 0.04958 1.018 0.309
Residual standard error: 18.18 on 1458 degrees of freedom
Multiple R-squared: 0.0007109, Adjusted R-squared: 2.55e-05
F-statistic: 1.037 on 1 and 1458 DF, p-value: 0.3086summary(m_year_trend)$r.squared[1] 0.0007108822m_year_trend60 <- lm(lay_DOY ~ year, data = timing_weather_60)
summary(m_year_trend60)
Call:
lm(formula = lay_DOY ~ year, data = timing_weather_60)
Residuals:
Min 1Q Median 3Q Max
-145.506 -10.576 -2.616 7.142 142.283
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 51.56562 99.65756 0.517 0.605
year 0.05049 0.04958 1.018 0.309
Residual standard error: 18.18 on 1458 degrees of freedom
Multiple R-squared: 0.0007109, Adjusted R-squared: 2.55e-05
F-statistic: 1.037 on 1 and 1458 DF, p-value: 0.3086Figures and Tables
cand_timing_30 <- cand_timing_30[names(cand_timing_30) != "null"]
cand_timing_60 <- cand_timing_60[names(cand_timing_60) != "null"]
build_obj3_base <- function(cand_list){
tab <- tibble::tibble(
model_id = names(cand_list),
Deviance = sapply(cand_list, stats::deviance),
K = sapply(cand_list, function(m) length(stats::coef(m))),
AICc = sapply(cand_list, AICcmodavg::AICc),
R2 = sapply(cand_list, function(m) summary(m)$r.squared),
`slope (SE)` = purrr::imap_chr(cand_list, ~{
broom::tidy(.x) |>
dplyr::filter(term != "(Intercept)") |>
dplyr::mutate(txt = paste0(
sprintf("%.3f", estimate),
" (", sprintf("%.3f", std.error), ")"
)) |>
dplyr::pull(txt) |>
paste(collapse = "; ")
})
) |>
dplyr::arrange(AICc) |>
dplyr::mutate(
`ΔAICc` = AICc - min(AICc),
w = exp(-0.5 * `ΔAICc`) / sum(exp(-0.5 * `ΔAICc`)),
Deviance = round(Deviance, 3),
AICc = round(AICc, 2),
`ΔAICc` = round(`ΔAICc`, 2),
w = round(w, 3),
R2 = round(R2, 3)
) |>
dplyr::select(
model_id, Deviance, K, AICc, `ΔAICc`, w, R2, `slope (SE)`
)
tab
}
model_key <- tibble::tribble(
~model_id, ~Group, ~Model,
"meanT_30", "Temperature", "Mean temperature (30 days)",
"maxT_30", "Temperature", "Maximum temperature (30 days)",
"minT_30", "Temperature", "Minimum temperature (30 days)",
"rainDays_30", "Rain", "Rainy days (30 days)",
"totalR_30", "Rain", "Total rainfall (30 days)",
"meanW_30", "Wind", "Mean wind speed (30 days)",
"maxW_30", "Wind", "Maximum wind speed (30 days)",
"meanT_60", "Temperature", "Mean temperature (60 days)",
"maxT_60", "Temperature", "Maximum temperature (60 days)",
"minT_60", "Temperature", "Minimum temperature (60 days)",
"rainDays_60", "Rain", "Rainy days (60 days)",
"totalR_60", "Rain", "Total rainfall (60 days)",
"meanW_60", "Wind", "Mean wind speed (60 days)",
"maxW_60", "Wind", "Maximum wind speed (60 days)",
"lagTotR", "Lagged effects", "Previous year's total rainfall",
"lagRainD", "Lagged effects", "Previous year's rainy days",
"lagMeanT", "Lagged effects", "Previous year's mean temperature",
"temp_rain", "Joint effects", "Mean temperature + rainy days",
"temp_wind", "Joint effects", "Mean temperature + mean wind speed",
"rain_wind", "Joint effects", "Rainy days + mean wind speed",
"temp_x_rain", "Interactive effects","Mean temperature × rainy days",
"minT_x_wind", "Interactive effects","Minimum temperature × mean wind speed",
"lagR_x_temp", "Interactive effects","Previous year's rainfall × mean temperature"
)
make_grouped_table <- function(cand_list, caption_text){
group_order <- c(
"Rain", "Temperature", "Wind",
"Lagged effects", "Joint effects", "Interactive effects", "Other"
)
tab <- build_obj3_base(cand_list) |>
dplyr::left_join(model_key, by = "model_id") |>
dplyr::mutate(
Group = dplyr::coalesce(Group, "Other"),
Model = dplyr::coalesce(Model, model_id),
Group = factor(Group, levels = group_order)
) |>
dplyr::arrange(Group, AICc) |>
dplyr::select(Group, Model, Deviance, K, AICc, `ΔAICc`, w, R2, `slope (SE)`)
group_sizes <- tab |>
dplyr::count(Group, name = "n") |>
dplyr::filter(n > 0) |>
dplyr::mutate(
start = cumsum(dplyr::lag(n, default = 0)) + 1,
end = cumsum(n)
)
kb <- knitr::kable(
tab |> dplyr::select(-Group),
caption = caption_text,
align = c("l","r","r","r","r","r","r","r","l"),
Twbooktabs = TRUE
) |>
kableExtra::kable_styling(full_width = TRUE)
for(i in seq_len(nrow(group_sizes))){
kb <- kb |>
kableExtra::pack_rows(
as.character(group_sizes$Group[i]),
group_sizes$start[i],
group_sizes$end[i],
bold = TRUE
)
}
kb
}
make_grouped_table(
cand_timing_30,
"Table 4. Models explaining laying date (day of year) at Dronfield using 30-day pre-laying weather cues. K = number of parameters; AICc = small-sample AIC; ΔAICc = difference from the best model; w = Akaike weight. R² is shown for all models; slope (SE) is reported for each weather predictor."
)| Model | Deviance | K | AICc | ΔAICc | w | R2 | slope (SE) |
|---|---|---|---|---|---|---|---|
| Rain | |||||||
| Rainy days (30 days) | 404580.3 | 2 | 12360.96 | 1024.96 | 0.000 | 0.161 | -2.103 (0.126) |
| Total rainfall (30 days) | 415698.5 | 2 | 12400.54 | 1064.54 | 0.000 | 0.138 | -0.286 (0.019) |
| Temperature | |||||||
| Minimum temperature (30 days) | 230069.4 | 2 | 11529.94 | 193.94 | 0.000 | 0.499 | -5.010 (0.131) |
| Mean temperature (30 days) | 241441.3 | 2 | 11600.33 | 264.33 | 0.000 | 0.475 | -5.478 (0.151) |
| Maximum temperature (30 days) | 280751.2 | 2 | 11820.41 | 484.41 | 0.000 | 0.389 | -4.775 (0.157) |
| Wind | |||||||
| Maximum wind speed (30 days) | 449433.4 | 2 | 12506.89 | 1170.89 | 0.000 | 0.022 | 4.154 (0.727) |
| Mean wind speed (30 days) | 452345.2 | 2 | 12516.32 | 1180.31 | 0.000 | 0.016 | 4.397 (0.915) |
| Lagged effects | |||||||
| Previous year's total rainfall | 444303.8 | 2 | 12278.29 | 942.29 | 0.000 | 0.010 | 0.013 (0.003) |
| Previous year's rainy days | 447834.4 | 2 | 12289.62 | 953.61 | 0.000 | 0.003 | 0.062 (0.032) |
| Previous year's mean temperature | 448801.1 | 2 | 12292.70 | 956.70 | 0.000 | 0.000 | -0.456 (0.576) |
| Joint effects | |||||||
| Mean temperature + rainy days | 203046.6 | 3 | 11349.65 | 13.65 | 0.001 | 0.558 | -5.016 (0.141); -1.508 (0.091) |
| Mean temperature + mean wind speed | 236458.4 | 3 | 11571.91 | 235.91 | 0.000 | 0.485 | -5.453 (0.150); 3.669 (0.662) |
| Rainy days + mean wind speed | 369830.8 | 3 | 12224.48 | 888.48 | 0.000 | 0.195 | -2.169 (0.120); 4.930 (0.828) |
| Interactive effects | |||||||
| Mean temperature × rainy days | 200878.1 | 4 | 11336.00 | 0.00 | 0.999 | 0.563 | -4.428 (0.204); 0.693 (0.563); -0.159 (0.040) |
| Previous year's rainfall × mean temperature | 235450.7 | 4 | 11373.62 | 37.62 | 0.000 | 0.476 | 0.026 (0.015); -4.875 (0.495); -0.001 (0.001) |
| Minimum temperature × mean wind speed | 215499.8 | 4 | 11438.51 | 102.51 | 0.000 | 0.531 | -1.416 (0.716); 11.636 (1.371); -1.198 (0.232) |
make_grouped_table(
cand_timing_60,
"Table 5. Models explaining laying date (day of year) at Dronfield using 60-day pre-laying weather cues. K = number of parameters; AICc = small-sample AIC; ΔAICc = difference from the best model; w = Akaike weight. R² is shown for all models; slope (SE) is reported for each weather predictor."
)| Model | Deviance | K | AICc | ΔAICc | w | R2 | slope (SE) |
|---|---|---|---|---|---|---|---|
| Rain | |||||||
| Rainy days (60 days) | 365166.2 | 2 | 12211.31 | 1450.61 | 0.000 | 0.243 | -1.350 (0.062) |
| Total rainfall (60 days) | 383366.2 | 2 | 12282.33 | 1521.62 | 0.000 | 0.205 | -0.160 (0.008) |
| Temperature | |||||||
| Minimum temperature (60 days) | 158946.3 | 2 | 10990.38 | 229.67 | 0.000 | 0.654 | -5.499 (0.105) |
| Mean temperature (60 days) | 172835.2 | 2 | 11112.60 | 351.89 | 0.000 | 0.624 | -6.105 (0.124) |
| Maximum temperature (60 days) | 209683.4 | 2 | 11394.57 | 633.86 | 0.000 | 0.544 | -5.834 (0.140) |
| Wind | |||||||
| Maximum wind speed (60 days) | 457772.5 | 2 | 12533.72 | 1773.01 | 0.000 | 0.004 | -1.889 (0.803) |
| Mean wind speed (60 days) | 459407.0 | 2 | 12538.92 | 1778.21 | 0.000 | 0.000 | -0.572 (0.998) |
| Lagged effects | |||||||
| Previous year's total rainfall | 444303.8 | 2 | 12278.29 | 1517.58 | 0.000 | 0.010 | 0.013 (0.003) |
| Previous year's rainy days | 447834.4 | 2 | 12289.62 | 1528.91 | 0.000 | 0.003 | 0.062 (0.032) |
| Previous year's mean temperature | 448801.1 | 2 | 12292.70 | 1531.99 | 0.000 | 0.000 | -0.456 (0.576) |
| Joint effects | |||||||
| Mean temperature + rainy days | 136023.4 | 3 | 10765.17 | 4.46 | 0.097 | 0.704 | -5.383 (0.116); -0.798 (0.040) |
| Mean temperature + mean wind speed | 171134.4 | 3 | 11100.18 | 339.48 | 0.000 | 0.628 | -6.150 (0.124); 2.328 (0.612) |
| Rainy days + mean wind speed | 336138.6 | 3 | 12085.12 | 1324.41 | 0.000 | 0.268 | -1.389 (0.060); -1.748 (0.855) |
| Interactive effects | |||||||
| Mean temperature × rainy days | 135421.3 | 4 | 10760.71 | 0.00 | 0.903 | 0.705 | -5.851 (0.217); -1.503 (0.280); 0.045 (0.018) |
| Previous year's rainfall × mean temperature | 167813.3 | 4 | 10889.01 | 128.31 | 0.000 | 0.626 | 0.067 (0.013); -4.329 (0.379); -0.004 (0.001) |
| Minimum temperature × mean wind speed | 153677.5 | 4 | 10945.22 | 184.51 | 0.000 | 0.666 | -2.223 (0.673); 11.131 (1.744); -1.092 (0.218) |
ggplot(timing_weather_60, aes(x = lay_DOY)) +
geom_histogram(bins = 15, fill = "grey70", colour = "black") +
labs(
x = "Laying date (day of year)",
y = "Frequency",
title = "Distribution of laying dates"
) +
theme_bw() +
theme(panel.grid = element_blank())
timing_weather_30_filt <- timing_weather_30 |>
dplyr::filter(year >= 1993, year <= 2024)
ggplot(timing_weather_30_filt, aes(x = year, y = lay_DOY)) +
geom_point(
size = 1.8,
colour = "grey30",
alpha = 0.6
) +
geom_smooth(
method = "lm",
se = TRUE,
colour = "grey",
linewidth = 1
) +
scale_x_continuous(
breaks = seq(1995, 2025, by = 5),
limits = c(1993, 2024)
) +
labs(
x = "Year",
y = "Laying date (day of year)",
title = "Annual variation in laying date at Dronfield"
) +
theme_bw() +
theme(
panel.grid = element_blank()
)+
geom_jitter(
width = 0.2,
height = 0,
size = 1.6,
colour = "grey30",
alpha = 0.6
)`geom_smooth()` using formula = 'y ~ x'Warning: Removed 45 rows containing missing values or values outside the scale range
(`geom_point()`).
# Figure 3.3: Temperature × rainfall interaction with prediction ribbons
m_int <- cand_timing_30$temp_x_rain
df <- timing_weather_30 |>
dplyr::filter(year >= 1993, year <= 2024)
# rainfall quantiles
rain_q <- quantile(df$rain_days_30,
probs = c(0.25, 0.5, 0.75),
na.rm = TRUE)
# prediction grid
newdat <- expand.grid(
mean_temp_30 = seq(
min(df$mean_temp_30, na.rm = TRUE),
max(df$mean_temp_30, na.rm = TRUE),
length.out = 80
),
rain_days_30 = as.numeric(rain_q)
)
# predictions + SE
pred <- predict(m_int, newdata = newdat, se.fit = TRUE)
newdat$fit <- pred$fit
newdat$lo <- pred$fit - 1.96 * pred$se.fit
newdat$hi <- pred$fit + 1.96 * pred$se.fit
newdat$Rain_group <- factor(
newdat$rain_days_30,
labels = c("Low rainfall (25%)",
"Medium rainfall (50%)",
"High rainfall (75%)")
)
# plot
ggplot(df, aes(x = mean_temp_30, y = lay_DOY)) +
geom_point(
colour = "grey75",
size = 1.6,
alpha = 0.35
) +
geom_ribbon(
data = newdat,
aes(
x = mean_temp_30,
ymin = lo,
ymax = hi,
group = Rain_group
),
alpha = 0.12,
inherit.aes = FALSE
) +
geom_line(
data = newdat,
aes(
x = mean_temp_30,
y = fit,
linetype = Rain_group
),
linewidth = 1.2,
inherit.aes = FALSE
) +
labs(
x = "Mean temperature (30 days pre-laying)",
y = "Laying date (day of year)",
linetype = "Rainfall conditions",
title = "Temperature–rainfall interaction predicts laying date"
) +
theme_bw() +
theme(
panel.grid = element_blank(),
legend.position = "right"
)
Laying date advanced with increasing pre-laying temperature, but this effect was contingent on rainfall conditions, indicating that thermal cues did not operate independently of broader environmental context. Warmer conditions are likely to reduce thermoregulatory costs and facilitate physiological readiness for breeding; however, the strength of this temperature effect varied with rainfall. Rather than rainfall acting as a direct trigger for breeding, wetter pre-laying conditions may reflect short-term improvements in ecosystem productivity or adult condition, thereby modifying how effectively temperature cues translate into earlier breeding.